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ABSTRACT 

We present a comparison between the observed color distribution, number and mass density of massive 
galaxies at 1.5 < z < 3 and a model by Hopkins et al. that relates the quasar and galaxy population 
on the basis of gas-rich mergers. In order to test the hypothesis that quiescent red galaxies are formed 
after a gas-rich merger involving quasar activity, we confront photometry of massive (M > 4x 10 10 Mq) 
galaxies extracted from the FIRES, GOODS-South, and MUSYC surveys, together spanning an area of 
496 arcmin 2 , with synthetic photometry from hydrodynamical merger simulations. As in the Hopkins 
et al. (2006b) model, we use the observed quasar luminosity function to estimate the merger rate. We 
find that the synthetic U — V and V — J colors of galaxies that had a quasar phase in their past match 
the colors of observed galaxies that are best characterized by a quiescent stellar population. At z ~ 2.6, 
the observed number and mass density of quiescent red galaxies with M>4x 10 10 Mq is consistent 
with the model in which every quiescent massive galaxy underwent a quasar phase in the past. At 
z ~ 1.9, 2.8 times less quiescent galaxies are observed than predicted by the model as descendants of 
higher redshift quasars. The merger model also predicts a large number and mass density of galaxies 
undergoing star formation driven by the merger. We find that the predicted number and mass density 
accounts for 30-50% of the observed massive star-forming galaxies. However, their colors do not match 
those of observed star-forming galaxies. In particular, the colors of dusty red galaxies (accounting for 
30-40% of the massive galaxy population) are not reproduced by the simulations. Several possible 
origins of this discrepancy are discussed. The observational constraints on the validity of the model 
are currently limited by cosmic variance and uncertainties in stellar population synthesis and radiative 
transfer. 

Subject headings: galaxies: evolution - galaxies: formation - galaxies: high redshift - galaxies: stellar 
content 
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1. INTRODUCTION 

In recent years, deep near- and mid-infrared observa- 
tions have revealed significant populations of red galaxies 
at redshifts z ~ 2 and above (Franx et al. 2003; Daddi 
et al. 2004; Yan et al. 2004). The population of Distant 
Red Galaxies (DRGs), selected by the simple observed 
color criterion J — K > 2.3, makes up 66% in number 
and 73% in mass of the 2 < z < 3 galaxy population at 
the high mass end (M > 10 11 Mq, van Dokkum et al. 
2006, see also Marchesini et al. 2007). Probing to lower 
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masses, Wuyts et al. (2007) found that the lower mass 
galaxies at redshifts 2 < z < 3.5 have bluer rest-frame 
U — V colors compared to the most massive galaxies. A 
substantial fraction of the massive red galaxies at high 
redshift are best characterized by a quiescent stellar pop- 
ulation on the basis of their broad-band SEDs (Labbe 
et al. 2005; Wuyts et al. 2007) and the presence of 
a Balmer/4000A break and absence of emission lines in 
their rest-frame optical spectra (Kriek et al. 2006). 

Any satisfying theory of galaxy formation has to ac- 
count for the presence and abundance of these mas- 
sive red galaxies in the early universe, a condition that 
was by no means met by the state-of-the-art hierarchical 
galaxy formation models at the time of their discovery 
(Somerville 2004). 

In the meantime, merger scenarios involving AGN 
activity have been invoked by semi-analytic models 
(Granato et al. 2004; Croton et al. 2006; Bower et 
al. 2006; De Lucia & Blaizot 2007; Somerville et al. 
2008) and hydrodynamical simulations (Springel et al. 
2005a; Di Matteo et al. 2005) to explain simultaneously 
the mass build-up of galaxies and the shutdown of star 
formation. Such an evolutionary scenario predicts an ob- 
scured (and thus red) star-burst phase and ends with a 
quiescent (and thus red) remnant galaxy (e.g., Hopkins 
et al. 2006a). Observational support for the connection 
between dust-enshrouded starbursts, merging, and AGN 
activity from samples of nearby Ultra-Luminous Infrared 
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Galaxies (ULIRGs) dates from as early as Sanders et 
al. (1988). Furthermore, the observed relation between 
the supermassive black hole (SMBH) mass and the mass 
(Magorrian et al. 1998) or the velocity dispersion (Fer- 
rarese & Merritt 2000; Gebhardt et al. 2000) of their host 
suggests that black hole and galaxy growth are intimately 
connected. This scaling relation can be reproduced by 
merger simulations with implemented AGN feedback (Di 
Matteo et al. 2005; Robertson et al. 2006b). 

Motivated by the observed and simulated correlations 
between the properties of SMBHs and their hosts, Hop- 
kins et al. (2006b) used the observed quasar luminosity 
function to derive the galaxy merger rate as a function 
of mass. This paper uses the merger rate function de- 
rived by this model in combination with hydrodynam- 
ical SPH simulations by Robertson et al. (2006a) and 
T. J. Cox to predict the color distribution, number and 
mass density of massive galaxies in the redshift range 
1.5 < z < 3 under the assumption that each galaxy once 
had or will undergo a quasar phase. We compare the re- 
sults to mass-limited samples in the same redshift inter- 
val, extracted from the multi-wavelength surveys FIRES 
(Franx et al. 2000; Labbe et al. 2003; Forster Schreiber 
et al. 2006a), GOODS-South (Giavalisco et al. 2004; 
Wuyts et al. 2008), and MUSYC (Quadri et al. 2007). 

The model we analyze in this paper resides in a much 
larger context that predicts that morphological transfor- 
mations, starbursts, quasars, and the growth of struc- 
ture are driven by galaxy mergers. This model has been 
well calibrated to many observations at low redshift (e.g., 
Jonsson et al. 2006 and Rocha et al. 2008 on attenuation 
of local spiral galaxies and mergers; Hopkins et al. 2008, 
2009 on the structure of local ellipticals) and high red- 
shift (e.g., Younger et al. 2009 and Narayanan et al. 2009 
on the infrared output of high-redshift (Ultra-)Luminous 
InfraRed Galaxies and Sub-Millimeter Galaxies). Here, 
we study the epoch of 1.5 < z < 3, when AGN and star 
formation activity was at its peak, and consider specifi- 
cally whether the observations of massive galaxies can be 
understood within this context. The comparison aims to 
shed light on the nature of massive galaxies and their evo- 
lutionary history, as well as identify where refinements to 
the model are needed. 

We give an overview of the observations and simula- 
tions in §2 and §3 respectively Next, the sample selec- 
tion is explained in §4. §5 addresses the methodology to 
populate a model universe with the binary merger simu- 
lations in order to predict number densities, mass densi- 
ties, and color distributions. We compare the predicted 
abundance of massive galaxies by the model to the ob- 
servations in §6. The optical and optical-to-NIR color 
distribution of observed and simulated massive galax- 
ies will be addressed in §7, followed by a discussion of 
their specific star formation rates (§8) and of the number 
and mass density of quiescent and star-forming massive 
galaxies in §9. We briefly compare observed and mod- 
eled pair statistics (§10) and address a few caveats on 
the observational and modeling results in §11. Finally, 
we summarize results in §12. 

We work in the AB magnitude system throughout this 
paper and adopt a HO — 70 km s~ x Mpc^ 1 , Q,m = 0.3, 
Oa = 0.7 cosmology. 

2. OVERVIEW OF THE OBSERVATIONS 



2.1. Fields, coverage, and depth 

We combine i^-band selected catalogs of three differ- 
ent surveys: FIRES, GOODS-South, and MUSYC. The 
reduction and photometry of the FIRES observations of 
the Hubble Deep Field South (HDFS) is presented by 
Labbe et al. (2003) and was later augmented with IRAC 
data. The field reaches a i^-band depth of 25.6 mag 
(AB, 5cr for point sources) and covers 5 arcmin 2 . It was 
exposed in the WFPC2 I/300, -8450, V^, Isu passbands, 
the ISAAC J s , H, and K s bands, and the 4 IRAC chan- 
nels. Following similar procedures, a i^ s -band selected 
catalog for the FIRES MS 1054-03 field was constructed 
by Forster Schreiber et al. (2006a). The field, covering 
19 arcmin 2 , has a -ft^-band depth of 25 mag (AB, 5<t for 
point sources). The catalog comprises FORS1 U, B, and 
V, WFPC2 Veoe, and I 81i , ISAAC J, H, and K s , and 
IRAC 3.6 /im - 8.0 pm photometry. 

Over a significantly larger area (114 arcmin 2 ), but to 
a shallower depth, a if s -band selected catalog, dubbed 
FIREWORKS, was constructed based on the publicly 
available GOODS-South data (Wuyts et al. 2008). The 
variations in exposure time and observing conditions be- 
tween the different ISAAC pointings lead to an inho- 
mogeneous depth over the whole GOODS-South field. 
The 90% completeness level in the i^-band mosaic is 
reached at an AB magnitude of K to t,AB = 23.7. The 
photometry was performed in an identical way to that 
of the FIRES fields, allowing a straightforward combi- 
nation of the three fields. The included passbands are 
the ACS -B435, V(;o6, *775, and zsso bands, the ISAAC 
J, H, and K s bands, and the 4 IRAC channels. We 
also use the ultradeep MIPS 24 /im (20 ^Jy, 5<r) imag- 
ing of the GOODS-South field. As for the IRAC bands, 
we used the information on position and extent of the 
sources from the higher resolution if s -band image to re- 
duce confusion effects on the 24 /im photometry (Labbe 
et al. in preparation). 

Finally, we complement the FIRES and GOODS-South 
imaging with optical-to-MIR observations of the MUSYC 
HDFS1, HDFS2, 1030, and 1255 fields for parts of our 
analysis. The X s -band selected catalogs, augmented 
with IRAC photometry, are presented by Marchesini et 
al. (2008). Together, the MUSYC fields span an area 
of 358 arcmin 2 . They reach the 90% completeness level 
at K to t,AB — 22.7. Given their shallower depth, they 
will only be used in the analysis of the most massive 
(M > 10 11 M Q ) high-redshift galaxies. 

2.2. Redshifts and rest-frame photometry 

Despite the large number of spectroscopic campaigns 
in the GOODS-South and FIRES fields, the fraction of 
-£T s -selected 1.5 < z < 3 galaxies that is spectroscopically 
confirmed is only 9%. The fraction drops to 3% when 
the MUSYC fields are included. Therefore, a reliable 
estimate of the photometric redshift is crucial in defining 
robust samples of massive high-redshift galaxies. 

Wuyts et al. (2008) used the EAZY photometric red- 
shift code by Brammer, van Dokkum & Coppi (2008) 
to fit a non-negative linear combination of galaxy tem- 
plates to the FIREWORKS J7 38 -to-8/zm spectral energy 
distributions of galaxies in the GOODS-South field. We 
applied an identical procedure to galaxies in the FIRES 
fields. The template set was constructed from a large 
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number of PEGASE models (Fiox & Rocca-Volmerange 
1997). It consists of 5 principal component templates 
that span the colors of galaxies in the semi-analytic 
model by De Lucia & Blaizot (2007), plus an additional 
template representing a young (50 Myr) and heavily ob- 
scured (Ay = 2.75) stellar population to account for the 
existence of dustier galaxies than present in the semian- 
alytic model. A template error function was applied to 
downweight the rest-frame UV and rest-frame NIR dur- 
ing the fitting procedure. The i\g-band magnitude m 
was used as a prior in constructing the redshift proba- 
bility distribution p(z\C, mo) for a galaxy with colors C. 
We adopt the value z mp of the redshift marginalized over 
the total probability distribution, 

_ J z p(z\C, to ) dz 
Zmp " Jp(z\C,m )dz ' W 

as best estimate of the galaxy's redshift. 

The uncertainties in the photometric redshifts were de- 
termined from Monte Carlo simulations. For each galaxy, 
a set of 100 mock SEDs was created by perturbing each 
flux point according to its formal error bar, and repeat- 
ing the Zphot computation. The lower and upper error 
on z P hot comprise the central 68% of the Monte Carlo 
distribution. 

We tested the quality of the photometric redshifts 
in two ways. First we compare them to the avail- 
able spectroscopic redshifts in the 1.5 < z < 3 inter- 
val, resulting in a normalized median absolute devia- 
tion (Jnmad ( Zp i+~^ ec ) = 0.075. The quality measure 

o~nmad remains the same when the spectroscopic red- 
shifts in the MUSYC fields arc included or excluded. 
Second we tested how well we could recover the red- 
shift from synthetic broad-band photometry of simulated 
SPH galaxies placed at redshifts 1.5 to 3. We found 
that the considered template set performed very well 
{o~nmad{Az/(1 + z)) = 0.027). The scatter in the com- 
parison to spectroscopically confirmed galaxies is larger 
than that derived from the simulations. This is likely 
due to the fact that the synthetic photometry is based 
on the same stellar population synthesis code as the tem- 
plate set used to recover the redshifts. Therefore, the 
second test only studies the impact of an unknown star 
formation history, dust and mctallicity distribution on 
the derived z p hot- 

We computed the rest-frame photometry by interpolat- 
ing between observed bands using the best-fit templates 
as a guide. Uncertainties in the rest-frame colors were de- 
rived from the same Monte Carlo simulations mentioned 
above, and comprise both a contribution from photo- 
metric uncertainties and from z p hot uncertainties. For a 
detailed description, we refer the reader to Rudnick et al. 
(2003). We used an IDL implementation of the algorithm 
by Taylor et al. (2009) dubbed "InterRest" . 

2.3. Stellar masses 

The stellar masses of the observed galaxies in FIRES 
and GOODS-South were derived by Forster Schreiber et 
al. (in preparation) following the procedure described 
by Wuyts et al. (2007). The stellar masses of galaxies in 
MUSYC were derived with the same method by March- 
esini et al. (2008). Briefly, we fit BC03 templates to 



the optical-to-8 ^m SED with the HYPERZ stellar pop- 
ulation fitting code, version 1.1 (Bolzonclla et al. 2000). 
We allow the following star formation histories: a sin- 
gle stellar population (SSP) without dust, a constant 
star formation history (CSF) with dust, and an exponen- 
tially declining star formation history with an e-folding 
timescale of 300 Myr (T300) with dust. The allowed Ay 
values ranged from to 4 in step of 0.2, and the attenu- 
ation law applied was taken from Calzetti et al. (2000). 
We constrain the time since the onset of star formation to 
lie between 50 Myr and the age of the universe at the re- 
spective redshift. Finally, we scale from a Salpeter (1955) 
IMF with lower and upper mass cut-offs of 0.1 M Q and 
100 M to a pseudo-Kroupa IMF by dividing the stellar 
masses by a factor of 1.6 (Franx et al. 2008). Whereas 
we adopt the BC03 models for our default analysis, we 
also performed the SED modeling with templates from 
Maraston (2005, hereafter M05) and otherwise identical 
settings. We indicate the results based on M05 models 
in the plots, and comment on them where relevant. On 
average, estimated stellar masses are lower by a factor 
1.5 when M05 models are used. 

2.4. Star formation rates 

We derived estimates of the total (unobscured plus ob- 
scured) star formation rate of the observed galaxies by 
adding the UV and IR light, scaled by the calibrations 
for the local universe (Kennicutt 1998): 

SFR [M Q yr- 1 } = 1.74 x lO" 10 [L m + 3.3L 2800 ) /L© 

(2) 

where the rest-frame luminosity L2800 = vL lJ (2SQQA) 
was derived from the observed photometry with the al- 
gorithm by Rudnick et al. (2003). The total IR lumi- 
nosity Ljr = L(8 — 1000 /xm) was derived from the ob- 
served 24 /xm flux density in combination with the pho- 
tometric redshift estimate (spectroscopic when available) 
following the prescription of Dale & Hclou (2002). As 
best estimate, we adopt the mean of the logarithm of 
all conversion factors corresponding to the Dale & Hclou 
(2002) IR spectral energy distributions within the range 
a. = 1 — 2.5, where a parameterizes the heating intensity 
level from active (a = 1) to quiescent (a = 2.5) galax- 
ies 1 . The variation from £jij ia= 2.5 to Lj_r. q= i is 0.9 dex 
in the redshift interval 1.5 < z < 3. Where relevant, 
we indicate this systematic uncertainty in the conversion 
from 24 /im to Lm and eventually star formation rate in 
the plots. 

3. OVERVIEW OF THE SIMULATIONS 

We use a set of smoothed particle hydrodynamics 
(SPH, Lucy 1977; Gingold & Monaghan 1977) simu- 
lations performed by Robertson ct al. (2006a) and 
T. J. Cox of co-planar and tilted, equal-mass, gas-rich 
(/gas = 0.8 at the start of the simulation) mergers over a 
range of galaxy masses. In §11, the validity of an equal- 
mass merger assumption is further discussed in the light 
of alternative mechanisms such as minor mergers (§11.5) 
and smooth accretion flows (§11.6). A description of the 
GADGET-2 code used to run the simulations is given 

1 We release the 24^tm-to-SFR conversion in table format on 
http: / / www.strw.leidenuniv.nl/fircworks 
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by Springel (2005b). Springel & Hernquist (2003) de- 
scribe the prescriptions for star formation and super- 
nova feedback. The interplay between the supermas- 
sive black hole(s) and the environment is discussed by 
Springel et al. (2005b). We refer the reader to Robert- 
son et al. (2006a) for specifications on this particular set 
of simulations and an explanation of how the progenitors 
were scaled to approximate the structure of disk galax- 
ies at redshift z = 3. The mass resolution varied from 
log mj ~ 5 per stellar particle for the lowest mass runs 
to log rrii ~ 6.5 per stellar particle for the most massive 
mergers. The photometry of the snapshots, stored with 
a time resolution of 14 Myr for the tilted and 70 Myr 
for the co-planar runs, was derived in post-processing as 
described by Wuyts et al. (2009). 

Briefly, the total attenuated spectral energy distribu- 
tion (SED) for a given snapshot consisting of N stellar 
particles is computed as follows: 



N 

£Att,tot(A) = m i ' Lint(agei, Z u A) 

»=1 



• exp 



a(X) 



(3) 



where m i; age^, and Zj are, respectively, the mass, age, 
and metallicity of stellar particle i that is treated as a 
single stellar population. Li nt is the intrinsic (unatten- 
uated) SED interpolated from a grid of templates from 
a stellar population synthesis code. Here, we use SSP 
templates from BC03 as default. Results obtained when 
using a grid of Maraston (2005, hereafter M05) SSP 
templates for different ages and metallicities will be ad- 
dressed as well. For the intrinsic emission Li nt of the 
black hole particle (s), we scale a luminosity-dependent 
template SED by the bolometric black hole luminosity 
given by the simulation (see Hopkins, Richards & Hern- 
quist 2007). Parameters in Eq. 3 that are dependent 
on the line of sight are subscripted with "los" . To each 
stellar particle, the column density of hydrogen (iVH Moe ) 
and the average metallicity along the line of sight (Z iy \ os ) 
was computed for 100 viewing angles, uniformly spaced 
on a sphere. The optical depth is proportional to this 
metallicity-scaled column density, so that A B /N m — 
(Z/0.02)(Ab/Nhi)mw where the gas-to-dust ratio of the 
Milky Way equals {A B /N m ) MW = 8.47 x 10~ 22 cm 2 . 
The wavelength dependence is adopted from an attenu- 
ation law (parameterized by the cross section c(A)). We 
use the Calzetti et al. (2000) reddening curve unless men- 
tioned otherwise. The change in predicted colors when 
adopting the SMC-like attenuation law from Pei (1992) 
will be discussed as well. 
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Fig. 1. — The relation between stellar mass and observed 
total K s magnitude for galaxies in the FIRES and GOODS- 
South fields at (a) 1.5 < z < 2.25 and (b) 2.25 < z < 3. The 
solid lines show the adopted log M > 10.6 (FIRES+GOODS- 
South) and logM > 11 (FIRES+GOODS-South+MUSYC) 
mass limits. The dotted line indicates the photometric limit 
of the GOODS-South imaging. The dashed line indicates 
the approximate limit for the MUSYC fields. There are few 
galaxies with logM > 10.6 and K 3t tot > 23.7, or logM > 11 
and -Ks, tot > 22.7. The largest incompleteness correction 
is needed for the highest redshift bin in the MUSYC fields. 
A fifth of the logM > 11 galaxies would be undetected 
by MUSYC, as estimated from the deeper FIRES+GOODS 
fields. 



4. SAMPLE SELECTION 

Our aim is to compare the color distribution, number 
and mass density of mass-limited samples of observed 
and simulated galaxies. We choose the mass-limit such 
that the observed sample is reasonably complete in the 
considered redshift interval, even for the field with the 
shallowest i\T s -band depth from which the sample was 
drawn. In order to optimally exploit the range in area 
and depth of the considered surveys, we define two mass- 
limited samples and divide each in two redshift bins: 
1.5 < z < 2.25 and 2.25 < z < 3, probing a similar 
comoving volume. The first sample contains galaxies 



more massive than logM = 10.6 (M ~ 4 x 10 10 M ) 
in the FIRES and GOODS-South fields. It contains 134 
and 106 objects in the low- and high-rcdshift bin respec- 
tively. We present the sample in Fig. 1, where we plot 
the stellar mass of all FIRES and FIREWORKS sources 
that are detected above the 5a level in the respective 
redshift bin against their total observed _ftT s -band magni- 
tude. The stellar mass correlates with the if s -band mag- 
nitude, but a scatter of an order of magnitude is present 
due to the range in redshifts and spectral types of the 
galaxies. The 90% completeness limit (K Sttot = 23.7) 
for the GOODS-South field , which is shallower than 
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the FIRES fields, is indicated with the dotted line. At 
1.5 < z < 2.25, no massive (logM > 10.6) galaxies 
fainter than K s x a t = 23.7 are found in the FIRES fields 
and deeper parts of the GOODS-South mosaic. The low- 
est i\~ s -band signal-to-noise ratio in the massive galaxy 
sample is S/Nk s — 9, strongly suggesting that no in- 
completeness correction is needed to compute the num- 
ber and mass density in the 1.5 < z < 2.25 rcdshift bin. 
In the 2.25 < z < 3 redshift bin, we find 7 well-detected 
massive (logM > 10.6) galaxies fainter than the 90% 
completeness limit of GOODS-South. All of them have 
5 < S/Nk s < 10, whereas the vast majority (90%) of 
massive galaxies with K s ^ ot < 23.7 are detected above 
the lOcr level. Evaluating the fraction of massive galaxies 
fainter than K Sitot = 23.7 in the area that is sufficiently 
deep to detect these sources, we estimate the complete- 
ness in the high-rcdshift bin to be ~ 93%. 

In order to reduce the uncertainty from cosmic vari- 
ance in the derived number and mass densities, we also 
compose a sample including the MUSYC fields, increas- 
ing the sampled area by roughly a factor of 3.6. The 
shallower depth forces us to restrict the mass limit to 
M > 10 11 Mq. The number of objects above this mass 
limit is 176 at 1.5 < z < 2.25 and 71 at 2.25 < z < 3.0. 
We derive the completeness in the two redshift inter- 
vals using the deeper FIRES and GOODS-South fields 
in Fig. 1. The dashed line marks the approximate depth 
(90% completeness) for the MUSYC fields. None of the 
1.5 < z < 2.25 galaxies with logM > 11 in the deeper 
FIRES and GOODS-South fields are fainter than this 
limit. For the 2.25 < z < 3 bin, the fraction of massive 
galaxies that would be missed by MUSYC increases to 
19%. In our analysis, we apply the appropriate incom- 
pleteness corrections unless stated otherwise. For each 
considered redshift bin and mass limit, the sample size 
decreases by roughly a factor 1.5 when using M05 mod- 
els. 

5. METHODOLOGY 

Large cosmological simulations with sufficient resolu- 
tion to study the accretion onto SMBHs are computa- 
tionally very challenging. First attempts were under- 
taken by Di Matteo et al. (2008), but by and large, 
hydrodynamical simulations including a self-consistent 
treatment of SMBH growth have only been run with ad- 
equate resolution on binary merger systems (Springel et 
al. 2005a; Di Matteo et al. 2005; Robertson et al. 2006a, 
2006b; Cox et al. 2006a, 2006b) or as zoom-in on over- 
dense regions of cosmological N-body simulations at very 
high redshift z ~ 6 (Li et al. 2007). In order to confront 
observations of 1.5 < z < 3 galaxies with the hydrody- 
namical simulations, we populate a model universe with 
the binary mergers according to a merger rate estimated 
from the observed quasar luminosity function following 
the prescription by Hopkins et al. (2006b). As a caveat, 
we note that this approach does not allow for a replen- 
ishment of the galaxy's gas reservoir by further accretion 
from the intergalactic medium (see also §11.5 and §11.7). 
As such, it is not a full cosmological prediction by itself, 
but our comparison can be used to see whether the as- 
sumption that the mergers will not experience further 
infall leads to reasonable results. 

Briefly, the conversion from quasar demographics to 
galaxy demographics goes as follows. From a large set 
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1 2 3 4 5 
z 

Fig. 2. — The birth rate of spheroids (in grayscale) as a func- 
tion of redshift and final stellar mass as derived from the observed 
quasar luminosity function. The meaning of the time scale arrow 
and the open circle is described in the text. The model by Hopkins 
et al. (2006b) assumes that this birth rate equals the merger rate of 
galaxies. The birth rate (i.e., merger rate) reaches a maximum of 
4.5 X 10~ 4 logM -1 Mpc~ 3 Gyr -1 at z ~ 2. As time evolves, the 
peak of the merger rate function shifts toward lower mass galaxies. 



of binary merger simulations, Hopkins et al. (2006a) de- 
termined the distribution of quasar lifetimes, describing 

the time rft ^^ p £^ fc ^ spent by a quasar of peak luminos- 
ity L pea k in the luminosity interval d\og(L). The ob- 
served quasar luminosity function simply corresponds to 
the convolution of this differential quasar lifetime with 
the birth rate n(L pea k) of quasars with peak luminosity 

Lpeak- 

®{L) = J dt ^g(^ fi( L P eak) d\ogL peak (4) 

Using a compilation of observed quasar luminosity 
functions in the hard X-rays (Ueda et al. 2003), soft 
X-rays (Hasinger, Miyaji,& Schmidt 2005), and optical 
(Richards et al. 2005), Eq. 4 was then de-convolved 
to solve for n(L pea k). The relation between peak lu- 
minosity of the quasar and the final black hole mass, 
derived from the same simulations, was then adopted to 
calculate the birth rate of black holes of a certain fi- 
nal mass 7i{Mbh)- This function was then converted 
to a birth rate of spheroids n{M sp h) as a function of 
their final stellar mass using the SMBH-host connec- 
tion M BH = 0.0012 , 5) 2 5 , M 8ph (Hopkins et al. 

( 1 + (TT75j ) 

2007). 

The model by Hopkins et al. (2006b) assumes that 
the birth rate of spheroids equals the major merger rate 
of galaxies. The resulting merger rate as a function of 
rcdshift and final stellar mass is displayed with grayscales 
in Fig. 2 (darker meaning a higher merger rate). Its 
rcdshift-dependence was derived by considering observed 
quasar luminosity functions at a range of redshifts. The 
peak of the merger rate at z ~ 2 has a value of 4.5 x 
10~ 4 logM -1 Mpc~ 3 Gyr^ 1 . A clear trend is visible 
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Fig. 3. — Typical evolution of a merger simulation: (a) star 
formation history, (b) history of the mass build-up, and (c) evolu- 
tion of the rest-frame V-band luminosity. The dotted line indicates 
when the peak in quasar luminosity is reached. For a detailed de- 
scription of the time evolution in these and other parameters (e.g., 
accretion rate, quasar luminosity, extinction) we refer the reader 
to Hopkins et al. (2006a). 



of mergers occurring in increasingly lower mass systems 
as we proceed in time (i.e., to lower redshifts) after this 
peak. If mergers are responsible for a significant part of 
the growth in stellar mass, this trend explains at least 
qualitatively the observed downsizing of star formation 
over cosmic time (Cowie et al. 1996). 

To evaluate the post-merger (i.e., post-quasar, since 
the merging event triggers quasar activity in the simu- 
lations) galaxy population at z ~ 2, we integrate the 
merger rate function from z = oo to 2 and over the 
whole stellar mass range. For example, when the in- 
tegration reaches (M*j ina i — 10 11 M Q ; z = 3), marked 
by the circle in Fig. 2, we compute the photometry of a 
merger simulation with a final stellar mass of 10 11 Mq 
at 1.1 Gyr after the peak of quasar luminosity (the time 
elapsed between z — 3 and z = 2). As explained in §3, 
we compute the synthetic photometry along 100 lines-of- 
sight, uniformly spaced on a sphere. The number density 
of galaxies at z = 2 with colors corresponding to the 100 
lines-of-sight is then scaled according to the value of the 
merger rate function at (M*ji na i = 10 11 M©; z = 3). 
Finally, a mass cut is applied to guarantee an identical 
selection of observed and simulated galaxies (§11.1.2 ad- 
dresses this step in more detail). 

In order to predict the abundance and properties of 
galaxies at z ~ 2 that have yet to reach their peak 
in quasar luminosity or did not even start merging at 
the evaluated epoch, one can in principle integrate the 
merger rate function down to lower and lower redshifts. 
How far one integrates beyond the evaluated redshift is a 
rather arbitrary choice. We caution that counting galax- 
ies long before they will contribute to the quasar lumi- 
nosity function will lead to large uncertainties given their 
unconstrained pre-merger history. The typical evolution 
of a merger simulation is illustrated in Fig. 3 where we 
plot the star formation rate, stellar mass, and rest-frame 



y-band luminosity as a function of time since the peak 
in quasar luminosity. We decide to integrate 700 Myr 
beyond the evaluated redshift, thus counting both the 
galaxies that are undergoing a merger-induced nuclear 
starburst (sometime between and 200 Myr before the 
quasar phase) and those with star formation triggered 
by the first passage (sometime between 200 and 700 Myr 
before the quasar phase). Hereafter, we will refer to all 
galaxies in an evolutionary stage between and 700 Myr 
before the quasar phase as merging galaxies. Such a pre- 
diction only counts those galaxies that will later merge 
and produce a quasar. Apart from predicting the abun- 
dance and properties of the post-quasar population, we 
will thus be able to constrain how many of the massive 
star-forming galaxies can be accounted for by merger- 
induced star formation. 

In the early stages of the merging process, the two 
progenitor galaxies may be resolved as two separate ob- 
jects in the observations. We simulate this in our model 
by calculating the projected angular distance between 
the central black holes for each snapshot and viewing 
angle at the considered redshift. If the projected sepa- 
ration is larger than l'/5, the 2 progenitors arc counted 
separately, each containing half the mass. For smaller 
projected separations, or when the black hole particles 
have merged, we consider the system unresolved and use 
the integrated properties in our model predictions. In 
cases where the two progenitor components are counted 
separately, the merger contributes twice to the number 
density, but with half the mass, and may therefore drop 
out of the mass-limited sample. 

Provided the assumption of a one-to-one correspon- 
dence between quasars and major mergers is valid, the 
formal uncertainty in the merger rate function presented 
in Fig. 2 originates mostly from the poorly constrained 
faint end of the observed quasar luminosity function, 
where one can assume a pure luminosity evolution or also 
a slope evolution. At the bright end, and therefore for 
our massive galaxy samples, the predictions are robust. 

6. THE NUMBER DENSITY, MASS DENSITY AND MASS 
FUNCTION OF GALAXIES WITH logM > 10.6 AT 
1.5 < Z < 3 

Before analyzing the observed and modeled massive 
galaxy sample as a function of color and galaxy type, 
we consider the overall abundance of galaxies above 
logM > 10.6. We computed the model number and 
mass density by integrating the merger rate function to 
700 Myr beyond the evaluated redshift, i.e., including 
galaxies up to 700 Myr before the quasar phase. In Fig- 
ure 4, the number and mass densities of galaxies with 
logM > 10.6 predicted by the model (empty symbols) 
are compared against the abundance of observed galax- 
ies (filled symbols) above the same mass limit. Circles 
indicate results using our default (BC03) SED modeling. 
Triangles represent the abundances derived using M05 
models. The results are listed in Table 1. The spread 
of the empty symbols indicates the freedom allowed by 
the model due to the poorly constrained faint end of the 
quasar luminosity function. 

We considered three sources of error in the observa- 
tions: Poisson shot noise, cosmic variance and selection 
uncertainties stemming from uncertainties in the redshift 
and mass estimates of individual galaxies. The black er- 
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Fig. 4. — The number and mass density of observed (filled sym- 
bols; FIRES + GOODS-S) and modeled (empty symbols) galaxies 
with logAf > 10.6 as a function of redshift. The black error bar 
represents the Poisson shot noise solely. The gray error bar ac- 
counts for uncertainties in redshift and mass, and a (dominating) 
contribution from cosmic variance. We find that both the pre- 
dicted number and mass densities agree within the error bars with 
the observed values. 



ror bars in Fig. 4 indicate the contribution from Poisson 
noise, ranging from 8 to 10%. We are more severely lim- 
ited by cosmic variance. We follow the method outlined 
by Somerville et al. (2004) to calculate the cosmic vari- 
ance as predicted from cold dark matter theory for a 
population with unknown clustering as a function of its 
number density and the probed comoving volume of the 
sample. The resulting contribution to the error budget is 
29% for the 1.5 < z < 2.25 and 30% for the 2.25 < z < 3 
redshift bin. Finally, the uncertainties in the individ- 
ual redshift and mass determinations propagate into the 
number and mass density of massive high-rcdshift galax- 
ies. We estimate the contribution to the total error bud- 
get from Monte Carlo simulations. We constructed 100 
mock catalogs for the FIRES and GOODS-South fields 
by perturbing the fluxes so that 68% of the perturbed 
values lie within the la errors. We then repeated the 
computation of photometric redshifts and other derived 
properties such as stellar mass. 

After constructing the 100 mock catalogs, we apply the 
same sample selection (redshift interval, log M > 10.6) 
and compute the number and mass density for each of 
them. The lower and upper limits comprising 68% of 
the distribution of mock number and mass densities were 
added in quadrature to the uncertainty from Poisson shot 
noise and cosmic variance, shown with the gray error 
bar in Fig. 4. The uncertainty in the number density 
propagating from redshift and mass uncertainties for in- 
dividual objects amounts to 5% and 10% for the low and 
high-redshift bin. The contribution to the uncertainty in 
the mass density is 6% and 14% for the low- and high- 
rcdshift bin respectively. We conclude that, even with 
the 138 arcmin 2 area of our combined deep fields, cosmic 
variance is still the limiting factor for the determination 
of the number and mass density of massive high-redshift 
galaxies. Furthermore, we note that the number and 




10.5 11.0 11.5 12.0 

log(M) [M ] 

Fig. 5. — The mass function of observed (red histogram; FIRES 
+ GOODS-S) and modeled (gray polygons) galaxies with log M > 
10.6 at redshift 1.5 < z < 2.25 (top panel) and 2.25 < z < 3.0 
(bottom panel). Merger remnants alone (t > *qso) cannot account 
for the total population of observed galaxies above the same mass 
limit. A better consistency with the observations is reached when 
integrating the merger rate function to include galaxies up to 700 
Myr before the quasar phase. 



mass densities systematically drop by a factor 1.5 and 
1.7 respectively when using M05 models. 

Given these uncertainties, Figure 4 shows a good agree- 
ment between the model number and mass density for 
the population of massive (logM > 10.6) galaxies as a 
whole and the observations. Plotting the mass function 
for the observations (red histogram) and the model (dark- 
gray polygon) in Figure 5, we find that the comparable 
abundance of observed and modeled galaxies still holds 
when studied as a function of galaxy mass. With lighter 
gray polygons, we illustrate the model prediction when 
including only galaxies up to 200 Myr before the merger 
(t > -200 Myr) or only merger remnants (t > tgso)- 
The width of the polygons reflects the uncertainty in the 
merger rate function. We conclude that merger remnants 
alone cannot account for the entire observed massive 
galaxy population at 1.5 < z < 3. However, including 
galaxies with merger-triggered star formation, the mass 
function predicted by the model is in good agreement 
with the observations. This encourages a more detailed 
investigation of the properties of observed and simulated 
massive galaxies. 

A detailed study of the stellar mass function of ob- 
served galaxies from z = 4.0 to z = 1.3, as well as a 
comprehensive analysis of random and systematic uncer- 
tainties, is presented by Marchesini et al. (2008). We 
note that our abundance estimates of observed massive 
galaxies are consistent with those derived from the work 
of these authors. Also, our cosmic variance estimates 
are consistent with those empirically derived by March- 
esini et al. (2008), who find that cosmic variance is the 
dominant source of random uncertainties at z < 2.5. 
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7. THE COLOR DISTRIBUTION OF GALAXIES WITH 
logM > 10.6 AT 1.5 < Z < 3 

7.1. The U — V color distribution 

First, we consider the optical color distribution of our 
sample of FIRES and FIREWORKS galaxies with M > 
4 x 10 10 Mq. A histogram of their rest-frame U — V 
colors is plotted in red in Fig. 6(a) and Fig. 6(b) for the 
low- and high-rcdshift bin respectively. No corrections 
for incompleteness were applied here, but we remind the 
reader that those are negligible for the low-redshift bin 
and 7% only for the high-redshift bin. 

Massive high-redshift galaxies span a broad U—V color 
range. In both rcdshift bins, the median color is U — V = 
1.5 and 68% of the galaxies lie within the 1.1 < U — V < 
1.8 interval. 

It is interesting to consider whether the simulated de- 
scendants and progenitors of quasars (or rather quasar 
hosts) above the same mass limit show colors that are 
similar and come in numbers comparable to those of the 
observed massive galaxy sample. In this section, we focus 
mainly on the first question, but note in passing that we 
show the predicted color distribution scaled to the same 
solid angle as probed by the FIRES and FIREWORKS 
observations. The filled gray histograms show the syn- 
thetic photometry of merger simulations in cither their 
post-quasar phase or in a phase of at most 700 Myr be- 
fore their peak in quasar luminosity. The numbers at 
each color are derived from the observed quasar lumi- 
nosity function by integrating the merger rate function 
from z = oo to 700 Myr beyond the evaluated redshift 
as described in §5. The colors of different evolutionary 
phases will be discussed separately in due time. The 
difference between the dark and light gray histogram re- 
flects the uncertainty in the merger rate function, itself 
due to uncertainties in the observed quasar luminosity 
function. Apart from uncertainties in the merger rate 
function, uncertainties in the synthetic photometry for 
a given simulation snapshot contribute to the total er- 
ror budget of the model predictions. To translate the 
simulated properties such as age, mass, and metallicity 
of the stellar particles to observables, we make use of 
a stellar population synthesis code to compute the in- 
trinsic colors and assume an attenuation law to calcu- 
late the dimming and reddening by dust. We investi- 
gate the dependence on attenuation law empirically by 
computing the synthetic photometry using a Calzetti et 
al. (2000) reddening curve and the SMC-like reddening 
curve from Pci (1992). We note that the synthetic colors 
derived with the Milky Way-like attenuation curve by 
Pei (1992) lie in between those produced by the two red- 
dening curves considered here. This is demonstrated by 
Wuyts et al. (2009). Similarly, we test the dependence 
on adopted stellar population synthesis templates empir- 
ically by computing the synthetic photometry based on a 
grid of BC03 single stellar populations (SSPs) and based 
on a grid of SSPs by M05. 

We note that the choice of attenuation law has a minor 
effect only on the U — V color. The use of M05 templates 
gives the simulated galaxies a slightly redder color. Over- 
all, the same conclusion can be drawn independent of the 
choice of attenuation law or stellar population synthesis 
code. Namely, the simulated galaxies with logM > 10.6 
span a color range that reaches from the bluest observed 



U-V colors to U - V ~ 2. At 1.5 < z < 2.25, the color 
distribution resembles remarkably well that of the bulk 
of the observed massive galaxies, both in shape and num- 
bers. We apply a Kolmogorov-Smirnov (K-S) test, and 
find that the observed and model color distributions do 
not differ at the 5% significance level. At 2.25 < z < 3, 
the predicted model colors do not reach the reddest U—V 
colors of observed galaxies above the same mass limit. A 
K-S test indicates a formal difference between the ob- 
served and model color distribution. We do note, how- 
ever, that the observed sources with U — V > 2 have 
fairly large uncertainties in their rest-frame optical color 
measurement, and are nearly all consistent within 1 a 
with an actual color of U — V < 2. The good overall 
correspondence between the observed and modeled op- 
tical color distributions gives a first indication that the 
number of massive post-quasar galaxies plus the number 
of galaxies in the process of merging at 1.5 < z < 3 as 
expected from the observed quasar luminosity function 
may account for a large fraction of the observed massive 
galaxy population in that redshift range. 

7.2. The V — J color distribution 

Turning to longer wavelengths, we now compare the 

V — J colors predicted for mergers and merger remnants 
(i.e., post-quasars) with masses above logM = 10.6 to 
the color distribution of observed galaxies in the same 
rcdshift interval and above the same mass limit (Fig. 7). 

Again, the color distribution of our observed massive 
galaxy sample has a large range of colors, reaching from 

V — J = 0.5 to V — J = 2.5. We measure a median 
color of V — J = 1.5 and 1.3 for the low- and high- 
rcdshift bin respectively. The central 68% intervals are 
1.1 < V - J < 1.9 and 0.9 < V - J < 1.8 

As for the U — V color distribution, we find that the 
adopted attenuation law has only a minor influence on 
the model color distribution, reaching at most shifts of 
0.2 mag toward redder V — J colors when the SMC-like 
reddening curve from Pei (1992) is used instead of the 
Calzetti et al. (2000) attenuation law. Comparing the 
model V— J color distribution derived from BC03 or M05 
templates immediately shows that the predictive power 
of the merger model is strongly hampered by the uncer- 
tainties in the rest-frame NIR wavelength regime that to- 
day's stellar population synthesis codes are facing. In the 
low- and high-redshift bin, the median V — J color of the 
model distribution is 0.4 and 0.5 mag redder when using 
M05 than when using BC03. One of the main differences 
between the BC03 and M05 templates is the treatment 
of thermally pulsating AGB stars. M05 uses the fuel 
consumption approach instead of the isochrone synthesis 
approach that BC03 follow. In addition, the two models 
construct the synthetic populations using different stellar 
isochrones. The combined effect is that M05 finds signif- 
icantly larger NIR luminosities for SSPs at ages between 
0.2 and 2 Gyr than BC03. For an in-depth discussion 
of the differences between the two codes, we refer the 
reader to Maraston (2005) and Maraston et al. (2006). 
It is worth stressing that, irrespective of whether the 
BC03 or M05 stellar population synthesis code is used, 
the red (V — J > 1.8) tail of the observed distribution 
has no counterparts in the modeled color distribution of 
merging and post-merger galaxies. Conversely, an excess 
of galaxies is found at blue (V — J <~ 0.9) or intermediate 
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Fig. 6. — The rest-frame U — V color distribution of observed galaxies with masses above logM = 10.6 in the FIRES and GOODS- 
South fields (solid and dashed line for masses based on BC03 and M05 respectively) for the redshift intervals (a) 1.5 < z < 2.25 and (b) 
2.25 < z < 3. Filled histograms are the predicted U — V color distribution of merging and post-quasar galaxies, scaled to the same solid 
angle as the observations. The light gray top of the model histogram reflects the uncertainty in the merger rate function. For a given 
redshift interval, the model predictions in the three panels give an indication of the uncertainty in the synthetic photometry induced by 
the choice of attenuation law (Calzetti et al. 2000 versus the SMC curve from Pei 1992) and the choice of stellar population synthesis 
code (BC03 versus M05). Overall, the predicted color distribution coincides with that of the observed massive galaxy sample, with roughly 
equal numbers. The red tail of the observed color distribution at 2.25 < z < 3 is not reproduced by the modeled merger and post-quasar 
population. 
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Fig. 7. — The rest-frame V — J color distribution of observed and simulated galaxies with masses above logM = 10.6 for the redshift 
intervals (a) 1.5 < z < 2.25 and (b) 2.25 < z < 3. Style as in Figure 6. The model V — J color distribution is only poorly constrained due 
to the uncertainties at NIR wavelengths in the stellar population synthesis codes. Nevertheless, we can conclude that there exist massive 
galaxies with redder V — J colors than those of modeled merging and post-quasar galaxies. 



(V - J ~ 1.4) optical-to-NIR colors for the BC03 and 
M05 model color distributions respectively. A K-S test 
shows that the difference between the distributions is sig- 
nificant at the 99.99% level. Uncertainties in the derived 
rest-frame colors of observed galaxies are unlikely to be 
responsible for the offset. 

7.3. U — V versus V — J color-color distribution 

7.3.1. Quiescent red galaxies 

Recently, a diagnostic color-color diagram of observer- 
frame I — K versus K-[A.b /im] has been proposed by 



Labbe et al. (2005) to distinguish three basic types of 
z > 2 galaxies. The rest-frame equivalent of this dia- 
gram, U — V versus V — J (hereafter UVJ) allows a com- 
parison of galaxies over a wider redshift range (Wuyts et 
al. 2007; Williams et al. 2009; Labbe et al. in prep). 
First, there are galaxies with relatively unobscured star 
formation, such as the majority of Lyman break galax- 
ies (Steidel et al. 2003) and their lower redshift BX/BM 
analogs (Adelberger et al. 2004). They typically have 
young ages and low reddening values, resulting in blue 
colors, both in the rest-frame optical and in the rest- 
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Fig. 8. — Model rest-frame U — V versus V — J color-color 
distribution of simulated galaxies with log M > 10.6 that have 
had a merger and quasar phase in their past {grayscales) , with 
a darker intensity indicating a larger number of post-quasars. 
Observed galaxies above the same mass limit in the FIRES 
and GOODS-South fields are overplotted. Red symbols mark 
the galaxies that satisfy the quiescent galaxy criterion, the 
selection window of which is marked by the red wedge. A 
notable difference between the synthetic photometry derived 
using the BC03 and M05 stellar population synthesis code is 
the redder V — J color in the latter case. Recognizing this 
uncertainty in the model prediction, we can still conclude 
that the predicted color distribution of post-quasars roughly 
coincides with that of quiescent red galaxies. 

frame optical-to-NIR. Second, there is a population of 
star-forming galaxies with much redder colors, due to 
the presence of dust. Their intrinsic (unobscured) colors 
are similar to those of Lyman break galaxies, but they 
are driven toward redder U — V and redder V — J colors 
along a dust vector whose slope depends on the nature 
and distribution of the dust (see, e.g., Wuyts et al. 2009). 
Finally, a population of galaxies with red U — V colors is 
present at z ~ 2 whose SED is well matched by that of 
a passive or quiescently star-forming galaxy at an older 
age. Their V — J colors are relatively blue compared to 
those of dusty starbursts at the same optical color. 

Labbe et al. (in preparation) designed a color criterion 
to select the quiescent red galaxies based on their rest- 
frame U , V, and J photometry. The selection window is 
defined as follows: 

U-V > 1.3 & V-J < 1.6 & U-V > 0.9(V-J)+X (5) 

where X is 0.31 and 0.4 for our low and high redshift 
bin respectively. 

The validity of this selection criterion was confirmed by 
the fact that galaxies within the UVJ box generally have 
low MlPS-based specific star formation rates (Labbe et 
al. in prep). Conversely, MIPS detected galaxies at z ~ 
2, suggesting dust-enshrouded star formation, tend to 
lie redward of the wedge. We draw the wedge in Fig. 
8 and indicate the location of all galaxies with logM > 
10.6 in the FIRES and GOODS-South fields in the color- 



color diagram. Red circles mark the objects that satisfy 
Eq. 5. The upper panels show the observed sample with 
BC03-based masses above logM > 10.6, whereas in the 
bottom panels the mass limit was applied to the M05- 
based stellar mass estimates. 

We also present a binned representation of the model 
color-color distribution of post-quasar galaxies only in 
Fig. 8. The panels correspond to the 1.5 < z < 2.25 
and 2.25 < z < 3 redshift bins, and model photometry 
derived from BC03 and M05 templates respectively. The 
color-color distribution computed with the SMC-like red- 
dening curve from Pei (1992) instead of the Calzetti et al. 
(2000) law is not plotted, but looks very similar. In the 
age range of 0.3 - 1 Gyr a dust-free SSP template of M05 
has redder V — J colors than a model of BC03 due to a 
different implementation of the thermally pulsing asymp- 
totic giant branch (TP-AGB) phase (e.g., Maraston et al. 
2006). Summing over all stellar particles (each treated 
as a SSP) of the simulated post-starburst (post-quasar) 
galaxies, this produces the shift toward redder integrated 
V — J colors in panels (c) and (d) with respect to panels 
(a) and (b). 

In all realizations of the synthetic photometry, the pre- 
dicted color-color distribution of the post-quasar popula- 
tion coincides more or less with the region of color-color 
space selected by the quiescent galaxy criterion. As a 
control sample, we analyzed a set of disk simulations, 
identical in initial conditions to the merger progenitors, 
but left to evolve in isolation instead of undergoing a vio- 
lent merger process. Due to the lack of new gas accretion, 
these simulated galaxies also reach phases of low specific 
star formation rate (SFR/M < 1/f Hubble)- However, 
the vast majority of these evolved disk galaxies have syn- 
thetic colors that place them outside the quiescent galaxy 
selection window, in a region of color-color space where 
only actively star-forming galaxies are found in the ob- 
servations. We conclude that an abrupt quenching of the 
star formation is required to reproduce the properties of 
observed quiescent galaxies at high redshift. Dissipative 
mergers involving AGN feedback seem to be a promising 
mechanism. 

7.3.2. Star-forming galaxies 

A significant fraction (~ two thirds) of the observed 
massive galaxy population at 1.5 < z < 3 has colors 
located outside the quiescent red galaxy wedge. These 
objects reach from blue U — V colors typical for Lyman 
break galaxies, which are known to host relatively un- 
obscured star formation, up to the redder optical and 
optical-to-NIR colors from galaxies that are believed to 
host heavily obscured star formation. Here, we inves- 
tigate whether the predicted color-color distribution for 
merging galaxies that will undergo a quasar phase in less 
than 700 Myr can reproduce the color range of observed 
star-forming galaxies. Fig. 9 compares the model predic- 
tion (grayscales) to the observed massive galaxy colors 
(blue circles for star- forming galaxies). 

As could be anticipated from §7.2, the model photome- 
try does not reproduce the colors of observed dusty star- 
forming galaxies (U — V > 1.3 and outside the quiescent 
red galaxy wedge). We note that our control sample of 
disk galaxies evolving in isolation do not reach the red 
optical-to-NIR colors of observed dusty starbursts either 
during their actively star-forming lifetime. 
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Fig. 9. — Model rest-frame U — V versus V — J color-color 
distribution of simulated galaxies with log M > 10.6 that will 
undergo a quasar phase in less than 700 Myr (grayscales), 
with a darker intensity indicating a larger number of galaxies. 
Observed galaxies above the same mass limit in the FIRES 
and GOODS-South fields are overplotted. Blue symbols mark 
the galaxies that fall outside the quiescent galaxy criterion 
(red wedge). A notable difference between the synthetic pho- 
tometry derived using the BC03 and M05 stellar population 
synthesis code is the redder V — J color in the latter case, (a) 
and (b) The model colors based on BC03 are a poor match to 
the observed star- forming galaxies (blue symbols). The V — J 
colors fall blueward of the observed distribution, and only the 
lower half of the observed U — V distribution of star-forming 
galaxies is reproduced, (c) and (d) The model colors based 
on M05 give a better match in the blue U — V regime, but 
simulated objects with V — J>1.6 are nearly absent. 



At bluer U — V, the synthetic photometry based on 
M05 templates gives a decent match to the observations, 
whereas the BC03 colors in combination with a Calzetti 
et al. (2000) attenuation law are offset by a few tenths 
of a magnitude toward bluer V — J. 

8. SPECIFIC STAR FORMATION RATE AS A FUNCTION 
OF STELLAR MASS 

So far, we have compared the synthetic colors of merg- 
ing and post-quasar galaxies with those of observed star- 
forming and quiescent galaxies respectively. The sepa- 
ration between star-forming and quiescent galaxies for 
our observed galaxies was based on their location in the 
UVJ diagnostic diagram. As an independent check, we 
now use the UV + 24 /im derived star formation rates 
to compare the observed distribution of specific star for- 
mation rates as a function of stellar mass with the distri- 
bution predicted by the merger model. The specific star 
formation rate is defined as the ratio of the star forma- 
tion rate over the stellar mass. It equals the inverse of 
a mass-doubling time in the case of constant star forma- 
tion. Here, we limit our sample to the GOODS-South 
field, where the 24 /im imaging is sufficiently deep (20 
/iJy, 5rj) to obtain useful constraints on the star forma- 
tion rates. 

Fig. 10 shows the binned model distribution in 



grayscales and overplotted are the observed massive 
galaxies that fall inside (red circles) and outside (blue 
circles) the quiescent red galaxy wedge. Upper limits are 
drawn for objects that were undetected by MIPS. Cross 
symbols mark those objects that are detected in the 1 
Ms Chandra X-ray exposure (Giacconi et al. 2002). We 
caution that the 24 /im flux of these objects could have 
an AGN contribution. Moreover, Daddi et al. (2007b) 
recently found that a significant fraction (20-30% to 
Kp* ga < 22, and up to - 50 - 60% for M ~ 10 11 M ) of 
star-forming galaxies that are not individually detected 
in the X-rays shows evidence for heavily obscured AGN 
by the presence of a mid-IR flux excess. The vertical er- 
ror bar indicates a conservative measure of the systematic 
uncertainty in the conversion from 24 /im flux to the ob- 
scured part of the star formation rate. The top and side 
panels show the distribution of masses and specific star 
formation rates separately. With lighter polygons, we 
illustrate how the predicted distribution changes when 
integrating the merger rate function only to the eval- 
uated redshift or 200 Myr past the evaluated redshift. 
The latter case includes the nuclear starburst phase, but 
not earlier star-forming phases. 

At 1.5 < z < 2.25, the broad-band color criterion is ef- 
ficient in distinguishing observed quiescent galaxies from 
star-forming galaxies with high specific star formation 
rates. In the higher redshift bin, we are more limited 
by upper limits on the 24 /im flux. The bulk of broad- 
band selected quiescent galaxies at 1.5 < z < 3 shows 
smaller specific star formation rates than their counter- 
parts outside the broad-band selection window, although 
some reach values above 1/tHubbic (dashed green line). 
The latter objects would double their stellar mass in less 
than a Hubble time if they form stars at a constant rate. 
Perhaps they are scattered into the UVJ box by pho- 
tometric redshift uncertainties, or their SFR is overesti- 
mated. 

The model SFR/M distribution is composed of 
merger-triggered star-forming galaxies with SFR/M > 
l/^Hubbio, and post-quasar systems with SFR/M < 
1/^Hubbic- The precise shape of the distribution at the 
low SFR/M end depends on, e.g., the relative orien- 
tation with which the disk progenitors were merged. 
However, the large difference between the star forma- 
tion mode before and after the quasar phase is a ro- 
bust feature of all considered simulations. The depth of 
the MIPS observations inhibits strong observational con- 
straints on the distribution of individual sources at the 
low SFR/M end, but a stacking analysis by Labbe et al. 
(in prep) reveals similarly low SFR/M values (~ 3 x 10~ 2 
Gyr -1 ) as for simulated post-quasar galaxies. As in the 
observations, in particular at 1.5 < z < 2.25, there is 
a slight hint that the most heavily star-forming objects 
reside primarily at the lower masses within our mass- 
limited sample. Papovich et al. (2006) and Reddy et 
al. (2006) find that the specific star formation rate is in- 
versely proportional to mass, implying that the ongoing 
star formation at z ~ 2 contributes more significantly to 
the mass buildup of low-mass galaxies than to high-mass 
galaxies. 

The predicted abundance of merger-triggered nuclear 
starbursts, occuring between and 200 Myr before the 
quasar phase, seems to be insufficient to account for all 
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Fig. 10. — Specific star formation rate as a function of stellar mass for massive galaxies at 1.5 < z < 3 in the GOODS-South field 
with colors falling inside (red circles) or outside (blue circles) the selection window for quiescent red galaxies. Cross symbols indicate 
which sources are detected in X-rays. The vertical error bar indicates the systematic error in SFR/M. The dashed green line marks 
1/tjjubblc- The model predictions are plotted with grayscales. The top and side panels show the mass and SFR/M distribution, with the 
black histogram representing the observed sample, and the grayscaled curves showing the model predictions for post-quasars and merging 
galaxies up to 700, 200, and Myr before the quasar phase. When integrating down to 700 Myr before the quasar phase, the predicted 
number density of galaxies with SFR/M > 1/tHubble i s 2-3 times smaller than observed, possibly (at least in part) due to AGN contribution 
to the 24 fim emission from which the observed SFR were derived. 



observed massive galaxies with high specific star for- 
mation rates (SFR/M > 1/tHubbie)- When including 
earlier phases of star formation induced by the merg- 
ing event (up to 700 Myr before the quasar phase), we 
find that the observed number density of galaxies with 
SFR/M > l/% u bbio is twice as large as predicted by 
the model. Part of this offset may be due to possible 
AGN contributions to the 24 fim emission from which 
the star formation rates were derived (see, e.g., Daddi et 
al. 2007b). Alternatively, this might indicate that not all 
star-forming galaxies can be accounted for by the consid- 
ered merger/quasar driven model for galaxy evolution. 

9. THE NUMBER AND MASS DENSITY OF MASSIVE 
GALAXIES AT 1.5 < Z < 3: ANALYSIS BY TYPE 

We now proceed to quantify the observed and modeled 
number and mass densities of different types of massive 
galaxies at 1.5 < z < 3. As before, the model prediction 
was derived by integrating the merger rate function to in- 
clude all galaxies that once contributed to the observed 
quasar luminosity function or will do so in less than 700 
Myr. From this, we extracted 6 samples using the cri- 
teria discussed in §7.3 and §8. Four are based on the 
UVJ diagnostic diagram: galaxies above log M > 10.6 
(or logM > 11) with broad-band colors satisfying the 
quiescent red galaxy criterion (Eq. 5, §9.1), and galaxies 
above logM > 10.6 (or logM > 11) that do not satisfy 
Eq. 5 (§9.2). Furthermore, we use the MlPS-based spe- 
cific star formation rates to independently address the 
abundance of relatively quiescent (SFR/M < 1/fHubbiei 
§9.1) and actively star-forming (SFR/M > 1/tHubbie, 
§9.2) galaxies above logM = 10.6. The logM > 11 sam- 
ples allow us to include the larger but shallower MUSYC 
survey in the comparison. In each case, we impose an 
identical selection criterion on the observed sample of 
galaxies. 



9.1. The number and mass density of massive quiescent 
red galaxies 

Having established the similarity in colors of the model 
post-quasar population and the observed quiescent red 
galaxy population above the same mass limit, we now 
turn to a comparison of their number and mass densi- 
ties. Our aim is to constrain the fraction (in number 
and mass) of massive quiescent red galaxies at redshifts 
1.5 < z < 3 that descendants of merger-triggered quasars 
can account for. 

In order to do this, we selected the observed and mod- 
eled galaxies with logM > 10.6 that lie inside the wedge 
defined by Eq. 5 and compute the number and mass den- 
sity for the probed comoving volume of ~ 3.5 x 10 5 Mpc 3 
in each redshift bin. The resulting number and mass den- 
sities are plotted as a function of central redshift of the 
redshift bin in Fig. 11(a). The filled circles and triangles 
indicate the observed number and mass density of qui- 
escent red galaxies above logM = 10.6 using BC03- and 
M05-based stellar masses respectively. Their values and 
corresponding uncertainties arc listed in Table 1. As in 
§6, the black error bars account for Poisson shot noise. 
The gray error bars also include selection uncertainties 
stemming from uncertainties in the redshift, mass, and 
rest-frame colors of individual galaxies, and a dominating 
contribution from cosmic variance. When correcting to 
the same IMF and mass limit, our estimates of the num- 
ber density are in good agreement with those of smaller 
samples of galaxies by Cimatti et al. (2008) and Kriek 
et al. (2008) that are spectroscopically confirmed to be 
quiescent. 

The empty symbols on Fig. 11(a) indicate the pre- 
dicted number and mass density of galaxies with log M > 
10.6 at 1.5 < z < 2.25 and 2.25 < z < 3 whose synthetic 
photometry places them within the selection wedge for 
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Fig. 11. — The number and mass density of observed (filled symbols) and modeled (empty symbols) massive galaxies as a function of 
rcdshift above the same mass limit and satisfying the same selection criterion. The black error bar represents the Poisson shot noise solely. 
The gray error bar accounts for uncertainties in redshift, mass, and rest-frame colors and a (mostly dominating) contribution from cosmic 
variance. The dashed error bars in panel (e) and (f) reflect the systematic uncertainty in the SFR. Arrows in panel (b), (d), and (f) 
indicate how the model results change when including the effect of mass underestimates in the SED modeling of star-forming galaxies. We 
consider 6 samples: (a) Quiescent red galaxies with logM > 10.6 in FIRES+GOODS-S, (b) Star-forming (non-quiescent) galaxies with 
logM > 10.6 in FIRES+GOODS-S, (c) quiescent red galaxies with logM > 11 in FIRES+GOODS-S+MUSYC, (d) star-forming galaxies 
with logM > 11 in FIRES+GOODS-S+MUSYC, and finally galaxies with (e) SFR/M < l/t Hub blo or (f) SFR/M > l/t Hu bble and 
logM > 10.6 in GOODS-S. We find that the model tends to overprcdict the number of quiescent and underpredict the number of actively 
star-forming galaxies by at most a factor 3. 
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Table 1. Number and mass densities for massive galaxies 



Observations' 1 Model Prediction' 



Type 


Mass limit 


Rcdshift 


n 




n 


P* 




M 




1(T 4 Mpc-' A 


10 7 M Q Mpc-' A 


1CT 4 Mpc-' A 


10 7 M Mpc' 3 


All 


4 x 10 10 


1.5 < z < 2.25 


4.0TM 


3.8±ii 


4.8 - 5.5 


5.4 - 6.0 


All 


4 x 10 10 


2.25 < z < 3 




9 O+0.9 
2 -°-0.9 


1.9 - 2.2 


1.9 - 2.2 


Quiescent 


4 x 10 10 


1.5 < z < 2.25 


i n+0.5 

^-o.b 


i 4+ - 6 


2.4 - 3.9 


3.0-4.8 


Quiescent 


4 x 10 10 


2.25 < z < 3 




n a+0-3 


0.4 - 1.1 


0.4 - 1.3 


Star-forming 


4 x 10 10 


1.5 < 2 < 2.25 


2.8ti; 


2 5+ 10 
z - d -0.9 


1.5 - 2.7 


1.2 - 2.6 


Star-forming 


4 x 10 10 


2.25 < z < 3 


9 O+0.9 


9 n +0.8 
z ' u -0.8 


1.0 - 1.7 


0.8 - 1.7 


Quiescent 


10 11 


1.5 < z < 2.25 


u -'-0.2 


1 1+0.3 
1 - 1 -0.3 


1.3 - 2.2 


2.3 - 3.6 


Quiescent 


10 11 


2.25 < z < 3 


0-2±g;i 


n 3 +0 - 2 

u -' 3 -0.1 


0.2 - 0.6 


0.3 - 1.0 


Star-forming 


10 11 


1.5 < z < 2.25 


o.8±g:i 


i 9+0.5 


0.3 - 0.9 


0.5 - 1.5 


Star-forming 


10 11 


2.25 < z < 3 


n 5+ - 3 

U.O_ 2 


o.8t8:I 


0.2 - 0.6 


0.3 - 1.0 


SFR/M < 1/tHubble 


4 x 10 10 


1.5 < z < 2.25 


i-3l° :S 


1 5+ - 7 


3.4 - 3.9 


4.2-4.8 


SFR/M < 1/tHubble 


4 x 10 10 


2.25 < z < 3 




n 5+ - 6 


1.0 - 1.2 


1.3 - 1.5 


SFR/M > 1/tHubble 


4 x 10 10 


1.5 < 2 < 2.25 


9 = + 1.0 
z - d -0.9 


9 9+O.9 
z - -0.8 


1.5 - 1.5 


1.1 - 1.2 


SFR/M > 1/tHubble 


4 x 10 10 


2.25 < z < 3 


2 6+ 1 ' 3 
z -°-1.0 




0.9 - 0.9 


0.7-0.7 



a The error bars in the observed densities account for Poisson noise, cosmic variance, and the uncertainties in rcdshift, rest-frame color and mass 
of the individual galaxies. They do not account for the systematic dependence on the stellar population synthesis code used to derive the stellar 
masses (values given here are for BC03), nor was the systematic uncertainty in the conversion from 24 fim to SFR (of the order of 1 dex) included 
in the results for the sample selected on SFR/M. 

^The range in model densities indicates a crude estimate of the size of uncertainties in the merger rate function and the dependence on choice of 
attenuation law and stellar population synthesis code to compute the synthetic photometry. 
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quiescent red galaxies. 95% of these modeled galaxies 
are in a post-quasar phase of their evolution. The differ- 
ent empty symbols represent predictions derived with the 
BC03 and M05 stellar population synthesis codes, with 
the Calzetti et al. (2000) attenuation law and the SMC- 
like attenuation law from Pei (1992). Their spread gives a 
crude indication of the freedom allowed by the model. It 
also takes into account the uncertainty in the merger rate 
function used to populate our model universe with the 
binary merger simulations. As noted already in §7.3.1, 
synthetic photometry of post-quasar galaxies based on 
M05 templates places part of them at redder V — J than 
the diagonal of the UVJ box, even though their sSFRs 
are low (10~ 2 — 10~ 05 Gyr -1 ). This explains why the 
M05-based model predictions in Fig. 11 (a) are lower than 
those based on our default BC03 models. We therefore 
consider the BC03-based model predictions of the quies- 
cent galaxy abundance as the most reliable. 

We find that the model predicts a number density of 
quiescent galaxies at z ~ 1.9 that is 2.5 times larger and 
a mass density that is 3 times larger than observed. At 
z ~ 2.6, the model and observations are consistent within 
the error bars. In other words, if anything, assuming 
a one-to-one correspondence between quasars and gas- 
rich mergers, the model by Hopkins et al. (2006b) over- 
predicts the abundance of merger remnants (i.e., post- 
quasar galaxies) at 1.5 < z < 3. The model predicts an 
increase by a factor 3.5 in the number and mass density 
for massive post-quasar galaxies in the 1 Gyr that passed 
between z = 2.6 and z — 1.9. The observed sample seems 
to suggest less evolution (a factor 1.2 in number density 
and 1.8 in mass density), and is even formally consistent 
with a non-evolving number and mass density over the 
considered redshift range. 

In order to reduce the effect of cosmic variance, we 
now include the MUSYC fields in our analysis, increas- 
ing the area by a factor 3.6 and reducing the cosmic vari- 
ance for a given mass limit with a similar factor. How- 
ever, this goes at the cost of depth: the 90% complete- 
ness limit for the MUSYC fields is 1 magnitude shal- 
lower than for GOODS-South. Consequently, we are 
restricted to a sample limited at M > 10 11 M©, even 
then requiring a 19% correction for incompleteness in 
the 2.25 < z < 3 bin. Since cosmic variance is larger 
for more massive galaxies, the reduction of uncertain- 
ties due to cosmic variance with respect to the deeper 
FIRES+GOODS-S sample is more modest than the in- 
crease in area. Fig. 11(c) shows the number and mass 
density of M > 10 11 M© galaxies that fall within the 
quiescent red galaxy wedge, as a function of redshift for 
the combined FIRES, FIREWORKS, and MUSYC sur- 
veys (filled black symbols) . Poisson noise is negligible for 
this sample. The gray error bars again account for cos- 
mic variance and the uncertainties in redshift, rest- frame 
color, and mass of the individual galaxies making up the 
sample. 

The main conclusion drawn from Fig. 11(c), based on 
a largely independent sample, is consistent with that of 
Fig. 11(a). Namely, a model in which every observed 
quasar produces a massive quiescent galaxy overpredicts 
the abundance of such galaxies at 1.5 < z < 3, by a 
factor of 3 for the sample above 10 11 M©. However, for 
our logAf > 11 sample we do find that the quiescent 
population grows by a similar amount (a factor <~ 3.5) in 



model and observations between z ~ 2.6 and z ~ 1.9. 

When comparing the abundance of quiescent galaxies 
selected by their MlPS-based SFR/M < l/i H ubbic (Fig. 
11(e)), we arrive at a similar conclusion. Here, we re- 
stricted ourselves to the GOODS-S field, for which deep 
MIPS 24/im data is available. If galaxies form stars at a 
constant rate, the adopted SFR/M threshold separates 
galaxies that double their mass in less than a Hubble 
time from those that formed the bulk of their stars prior 
to the time of observation. We note that this threshold 
coincides with the minimum in the observed and mod- 
eled SFR/M distribution (Fig. 10). For the high-redshift 
(z ~ 2.6) bin, the presence of a bimodal SFR/M distri- 
bution could not be observationally confirmed or ruled 
out due to the larger upper limits on individual SFR/M 
measurements. As in Fig. 11(a), we find the model num- 
ber and mass density to be larger than observed, by a 
factor 2.5 and 3 respectively. 

Possibly, our results suggest that part of the quasar- 
descendants are rejuvenated, e.g. by new gas infall fu- 
eling recurrent star formation. Large cosmological simu- 
lations including a prescription for AGN feedback, such 
as undertaken by Di Matteo et al. (2008), but with a 
resolution comparable to the binary merger simulations 
presented here are necessary to investigate such scenar- 
ios. 

We should also bear in mind that the evolution of 
the Mbh — M* relation and its scatter with redshift is 
only poorly constrained. Since this relation is used to 
translate observed quasar demographics into galaxy de- 
mographics (see §5), an uncertainty of a factor of a few in 
Mbh — Af* might resolve the offset in abundance as well. 
Finally, we caution that beaming effects at the bright 
end of the quasar luminosity function might complicate 
the methodology described in §5. 

9.2. The number and mass density of massive 
star-forming galaxies 

Following identical procedures as outlined above, we 
analyze the number and mass density of massive galaxies 
with colors outside the quiescent red galaxy wedge in 
Fig. 11(b). Again, we used a Monte Carlo simulation 
to determine how many galaxies moved into or out of 
the selection window when perturbing their photometry 
within the error bars and hence changing the derived 
properties such as mass and rest-frame photometry. 

We find massive star-forming galaxies in the observed 
fields to be 2.3 times more abundant in number and 2 
times more in mass than massive quiescent galaxies at 
1.5 < z < 3. Their contribution to the overall number 
and mass density decreases slightly, to 60%, when con- 
sidering our logM > 11 sample (Fig. 11(d)). 

Given that M05-based synthetic photometry placed 
some of the simulated post-quasar galaxies with low sS- 
FRs in the star-forming part of the UVJ diagram, we 
focus on the results obtained using BC03 templates (cir- 
cles). The model seems to predict that massive star- 
forming systems are less, rather than more, abundant 
than high-redshift quiescent systems above the same 
mass limit. Whereas the model overpredicted the quies- 
cent number and mass density by factors of a few (§9.1), 
the abundance of the star-forming population is under- 
predicted by a similar amount: a factor 1.8 at z ~ 1.9 
and a factor 2.2 at z ~ 2.6. 
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We now compare the abundance of actively star- 
forming galaxies again, using a sSFR threshold as se- 
lection criterion, rather than the rest-frame optical to 
NIR colors. Selecting galaxies with SFR/M > 1/^Hubbie 
from FIREWORKS, we find an observed number den- 
sity of 2.5 x 1CT 4 Mpc~ 3 and 2.6 x KT 4 Mpc~ 3 at 
1.5 < z < 2.25 and 2.25 < z < 3 respectively (Fig. 11(f)). 
Since we interpreted all the 24 fim emission as dust re- 
emission from star formation, the true number density 
can be lower depending on the contribution from AGN 
(see, e.g., Reddy et al. 2005; Papovich ct al. 2006; Daddi 
et al. 2007b). The merger model predicts an abundance 
of galaxies with high specific star formation rates that is 
1.7 (3.0) times smaller than observed in our low (high) 
rcdshift bin. Despite the possible AGN contribution to 
the 24 /jm emission and the large systematic uncertainty 
in the conversion from 24 /zm to the dust-obscured con- 
tribution to the star formation rate (dashed line in Fig. 
11(f)), the different selection of actively star-forming sys- 
tems produces a similar result as derived from the UVJ 
diagram. 

Testing SED modeling on mock observations of hydro- 
dynamic merger simulations, Wuyts et al. (2009) cau- 
tion for systematic mass underestimates during phases 
of merger-induced star formation. As discussed in §11.1, 
accounting for such an effect would increase the offset 
with respect to the observed abundance, such that the 
merger model then accounts for a third of the observed 
star-forming galaxies with logM > 10.6 at z ~ 1.9, and 
a quarter at z ~ 2.6. This result may be consistent with 
recent evidence for other mechanisms than major merg- 
ers driving part of the star formation at high rcdshift 
(Daddi et al. 2007; Shapiro et al. 2008; Genel et al. 
2008). We note that the fraction of star- forming galaxies 
accounted for by the merger model decreases when only 
considering the short-lived (~ 100 Myr) starburst at final 
coalescence instead of counting all merger-triggered star- 
forming phases as we do by integrating the merger rate 
function to 700 Myr before the peak of quasar activity. 

Overall, we conclude that the model abundances of 
massive galaxies are formally consistent with the ob- 
servations. However, dividing the sample in quiescent 
and star-forming populations, offsets of factors 2-3 occur, 
with the model predictions being higher for quiescent and 
lower for star-forming galaxies. 

10. PAIR STATISTICS 

In our model predictions of the number and mass den- 
sities of different samples of massive galaxies, a simple 
criterion determines whether we use the integrated prop- 
erties (color, mass, SFR) for the merging pair or treat 
the two progenitors as resolved systems, counting each 
one separately, contributing half the mass and SFR (see § 
5). Evolutionary phases, redshifts and viewing angles for 
which the projected angular separation between the cen- 
tral SMBHs is less than 175 were considered unresolved. 
All post-quasar predictions are robust against the precise 
form of the criterion, since by that time the two progeni- 
tors have formed one galaxy. For earlier phases, applying 
the criterion decreases the mass density of massive galax- 
ies, since galaxies drop out of the mass-limited sample. 
The effect on the number density is less trivial. On the 
one hand, galaxies drop out of the mass-limited sample. 
On the other hand, some merging pairs contribute twice. 
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Fig. 12. — Relative galaxy density (<5) as a function of massive 
(logM > 10.6) 1.5 < z < 3 galaxy to z > 1.3 galaxy separation 
(Rs) in the GOODS-South field. The distribution predicted by 
the merger model is indicated with the dashed line. At separations 
smaller than 1'.'5 (dotted line) an increasing number of galaxy pairs, 
if present, will be missed because they would be detected as a 
single object. We find a clear excess at small pair separations 
(R s < 8"), as predicted by the merger model. A weak pair excess 
is also visible when only considering the distribution of separations 
between massive 1.5 < z < 3 galaxies (inset panel), but the excess 
is much below the prediction. 



Here, we focus on an additional test of the merger 
model allowed by the fact that some of the pairs will 
be resolved into two objects. If a significant fraction of 
the massive galaxy population at 1.5 < z < 3 is indeed 
related to merging events, as our analysis suggests, we 
expect to see an excess in the pair statistics with respect 
to a random distribution of galaxies on the sky. 

We present the distribution of galaxy-galaxy separa- 
tions in the GOODS-South field in Figure 12 (solid his- 
togram). We decide not to include the other fields, to 
prevent differences in depth from influencing the pair ex- 
cess signal. The main panel shows the results from a 
cross-correlation of our massive (logM > 10.6) galaxy 
sample at 1.5 < z < 3 with the sample of all galaxies 
above z > 1.3 in the GOODS-South field, thus avoiding 
the risk of losing pair members that by a typical photo- 
metric redshift error were placed at some lower redshift. 
For each massive galaxy at 1.5 < z < 3, we measure the 
distance to all z > 1.3 galaxies. We compute the statistic 
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where j is the total number of objects in our massive 
galaxy sample and Ni(R s ) is the number of z > 1.3 galax- 
ies that lie between a distance R s — e and R s + e from 
galaxy i. For a random uniform distribution of galaxies, 
S(R S ) will be flat. Figure 12 shows that for our sample 
of massive galaxies at 1.5 < z < 3, this is clearly not 
the case. An excess of pairs at R s < 8" is visible, also 



when we consider the distribution of separations between 
members of the massive galaxy sample at 1.5 < z < 3 
only (inset panel). We note that, using the wider area 
UDS field, Quadri et al. (2008) found an upturn of the 
correlation function of distant red galaxies at 2 < z < 3 
on similarly small scales (9 < 10"). Spectroscopic rel- 
ative velocity measurements are needed to assess what 
fraction of the small scale excess is due to bound pairs in 
the process of merging, and how much can be attributed 
to an enhanced number of projected pairs due to cluster- 
ing. 

From the simulations, we measured the physical sep- 
arations between the 2 merging galaxies and computed 
the distribution of separation angles in arcseconds on the 
sky using the merger rate function (see §5) and assum- 
ing random viewing angles. Adding the mean value of S 
as measured in the interval 30" < R s < 80", we obtain 
a model prediction (dashed line) that is in qualitative 
agreement with the cross-correlation results, but larger 
than the weak pair excess seen in the auto-correlation. 
Admittedly, the predicted distribution is subject to the 
orbital configuration set at the start of the simulation, 
an effect that is not explored in this paper. 

11. COMMENTS AND CAVEATS 

In this Section, we list a number of caveats, and indi- 
cate prospects for improvements on both the model and 
observational side. We discuss aspects affecting the de- 
termination of galaxy abundances as well as a number 
of possible reasons for the discrepancy between the syn- 
thetic and observed colors of massive star-forming galax- 
ies. Future investigations along these lines will help to 
further test the merger model. 

11.1. Simulating the observing procedure 
11.1.1. Color gradients 

First, it is possible that the colors of observed and mod- 
eled galaxies are in fact in agreement, but that a discrep- 
ancy was found because we did not simulate the whole 
observing procedure. The observed optical-to-NIR colors 
are measured on PSF-matched images within apertures 
of size 1" to 2" . The IRAC photometry was performed 
within apertures of 3" diameter, and scaled to the smaller 
color aperture assuming all sources had a flat K s - IRAC 
color profile (see Wuyts et al. 2008). The synthetic colors 
instead were based on integrated photometry of all stel- 
lar particles, irrespective of their distance to the galaxy 
center. The presence of a color gradient with redder emis- 
sion in the central regions of the galaxies could therefore 
induce an offset in colors in the observed direction. 

In Fig. 13, we investigate the presence of color gra- 
dients for three snapshots (before, during and after the 
merger) of a simulation with final stellar mass M*fi na \ = 
1.2 x 10 11 Mq. The color measured within a radius r is 
plotted as a function of that radius. The polygon for a 
given snapshot illustrates the range in colors as we view 
the galaxy from different angles. The size of the aper- 
tures used for the observations is indicated for reference. 

A color gradient is clearly present in all three phases, 
with the color getting progressively bluer as we increase 
the aperture size. The gradient is most pronounced dur- 
ing the merger event, and more so for V — J than for 
U — V. During the nuclear starburst (green polygon), the 
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Fig. 13. — Simulated rest-frame U — V and V — J colors be- 
fore (blue), during (green), and after (red) a merger that ends with 
final — 1.2 X lO 1 ^ Mq as function of radius in which the color 
was measured. The height of the polygons indicates the range of 
colors for the same snapshot, as seen from different angles. Aper- 
ture sizes for the optical-to-NIR (1" < diameter < 2") and IRAC 
photometry (3" diameter) of real galaxies are marked in gray. Color 
gradients are present with the galaxy core being redder than the 
outskirts, particularly during the merger-driven nuclear starburst 
(green). 



difference (for a 1"5 diameter aperture) amounts up to 
0.15 mag in U — V. Given that our V — J estimates were 
effectively derived from 3" diameter apertures, a more so- 
phisticated simulation of the observing procedure could 
redden the V ~ J colors by up to 0.1 mag. For a given 
aperture size, the color difference between the aperture 
and integrated photometry also increases as we consider 
simulations of larger mass. E.g., for a merger that is five 
times more massive than that shown in Fig. 13, we find 
(U — V) d _ 1 n 5 - (U — V)totai reaches up to 0.25 mag and 
(V - J)d=3>> - (V - J)totai up to 0.6 mag. 

We note that, since color apertures varied from ob- 
ject to object in the observations, incorporating the de- 
tails of the observing procedure in our comparison is not 
straightforward. Moreover, it would require a proper 
treatment of smearing by the PSF, which was not ap- 
plied in Fig. 13. Nonetheless, our analysis indicates that 
color profiles are present in various degrees during differ- 
ent phases of the merger scenario, and can contribute 
to the photometric differences between simulated and 
real galaxies. A study of color profiles for high-redshift 
galaxies of different types, such as will be enabled by the 
high resolution imaging with WFC3 onboard HST, will 
provide additional constraints on the role of mergers in 
galaxy evolution. 

11.1.2. Accounting for biases in SED modeling 

In comparing abundances of modeled and observed 
galaxies, it is critical to apply identical selection crite- 
ria to both samples. In this paper, we work with mass- 
limited samples, where the stellar mass of modeled galax- 
ies is known from the simulation output, and that of 
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the observed galaxies is derived by fitting stellar popula- 
tion synthesis templates to their multi- wavelength SED. 
Wuyts et al. (2009) tested the performance of the latter 
method by computing mock high-redshift observations 
of the same simulations as used in this paper, feeding 
them to the same SED modeling procedure as described 
in §2.3, and comparing the derived masses with the true 
values known from the simulation output. The mass es- 
timates of merger remnants are very robust. However, 
during phases of merger-induced star formation system- 
atic mass underestimates of 0.1 - 0.2 dex occur, with a 
tail toward more severe underestimates. We repeated 
our comparison of galaxy abundances, not using the true 
mass of the modeled galaxies but the value derived from 
fitting their virtually observed SEDs. As expected from 
the results by Wuyts et al. (2009), the modeled abun- 
dance of massive quiescent galaxies remains nearly un- 
affected. A significant fraction of simulated star-forming 
galaxies drops out of the mass-limited sample, reducing 
the modeled number density by a factor 1.6. The mass 
density decreases by a factor 1.8. The shifts in number 
and mass density when including the effects of mass un- 
derestimates of star-forming galaxies in SED modeling 
are indicated with arrows in Fig. 11. We conclude that, 
accounting for biases in SED modeling, the merger model 
predicts that merger-induced star-forming galaxies above 
logM > 10.6 can account for a third (quarter) of the to- 
tal observed massive star- forming population at z ~ 1.9 
(2.6). Similar merger fractions among high-redshift star- 
forming galaxies are found by kinematic studies of the 
SINS survey (Shapiro et al. 2008). 

11.2. Dependence on stellar population synthesis and 

IMF 

As pointed out in §7.2, the predicted rest-frame NIR 
luminosities for a single stellar population of a given mass 
are brighter for the M05 than for the BC03 stellar popu- 
lation synthesis code. Consequently, the mass estimates 
for observed galaxies with ages between 0.2 and 2 Gyr 
are lower by about a factor 1.5 when modeled with M05 
instead of BC03 templates. We indicated the resulting 
systematic uncertainties in the observed and modeled 
number and mass densities. We find number densities 
for all samples discussed in this paper to be typically 
two thirds and mass densities to be ~60% of the value 
obtained with BC03 masses. Whereas using M05-based 
masses brings the observed abundance of star-forming 
galaxies in agreement with the model prediction, it in- 
creases the discrepancy between model and observations 
for the quiescent population. Furthermore, deviations 
from a Kroupa (2001) IMF would change our results 
on number and mass densities of massive galaxies. Re- 
cently, van Dokkum (2008) and Dave (2008) presented 
evidence for an evolving IMF that is more weighted to 
high mass stars at higher redshift. For the mass limit of 
log M > 10.6 considered in this paper, the general trend 
of such an evolving IMF would be to lower the abundance 
of high-redshift galaxies above a certain mass limit, but is 
not trivial to implement since the change in mass-to- light 
ratio would depend on the galaxy's age. Marchesini et al. 
(2008) arc the first to investigate the effect of bottom- 
light IMFs on the stellar mass function, and find that it 
does not merely result in a shift. The precise shape of the 
stellar mass function depends on the characteristic mass 



of the bottom-light IMF, and can in some cases, perhaps 
countcrintuitively, lead to abundances at the very high- 
mass end (logM > 11.5) that exceed those derived with 
a standard IMF. 

Concerning the discrepancy in colors, this could be due 
to an incorrect modeling of the stellar populations, rather 
than invalid assumptions at the basis of the model (i.e., 
the one-to-one correspondence between quasars and gas- 
rich mergers). Apart from the choice of stellar popu- 
lation synthesis code (see §7.2), the synthetic photom- 
etry depends on the attenuation law applied to each of 
the stellar particles. We note however that the use of a 
Milky Way-like attenuation law from Pei (1992) leads to 
colors intermediate between those based on the Calzetti 
et al. (2000) and SMC-like (Pei 1992) attenuation laws 
presented in this paper. An attenuation law that is less 
gray than that of the SMC would be required to repro- 
duce redder colors for dusty starburst galaxies. 

Another stellar population parameter influencing the 
synthetic colors is the metallicity of the gas and the stars. 
In this paper, we adopted initial gas metallicitics derived 
from the closed box model (Talbot & Arnett 1971) for 
the 80% gas fraction (/ gas ) at the start of the simulation: 

Zinit = -t/ln(/ gas ) (7) 

where y = 0.02 is the yield. The simulation keeps track of 
the subsequent evolution in the gas metallicity, and stel- 
lar metallicitics are based on the metallicity of the gas 
out of which they form. If the gas was pre-enriched, this 
would boost the optical depths and redden the colors. 
Evidence of high (~ Zq) metallicitics of massive high- 
redshift galaxies with red colors is given by van Dokkum 
et al. (2004). Repeating the post-processing of simu- 
lation snapshots with IZq added to the gas and stellar 
metallicities, we obtain colors that are 0.1 to 0.4 mag 
redder in U — V and 0.1 to 0.9 mag redder in V — J. We 
note however that in V — J the largest increase occurs for 
blue galaxies and the color distribution based on BC03 
does not reach beyond V — J ~ 1.8. 

11.3. Radiative transfer 

The fact that we find a difference between the observed 
and modeled V — J colors of star-forming galaxies, does 
not necessarily mean that the merger scenario (mergers 
triggering starbursts and quasar activity, and leaving a 
quiescent remnant) should be abandoned. Neither does 
it have to imply a failure of the hydro-simulations to 
realistically model such a merger event. Instead, and 
in fact more plausibly, it might reflect the difficulty of 
translating physical parameters stored in the GADGET- 
2 output into observables such as colors and fluxes. So 
far, we followed Hopkins et al. (2005a) and Wuyts et al. 
(2009) to compute the synthetic photometry assuming 
the cold gas clouds have a negligible volume filling fac- 
tor and ignoring the effect of scattering. In this subsec- 
tion, we explore how the synthetic photometry changes 
when adopting the recently developed radiative transfer 
code SUNRISE (Jonsson 2006). SUNRISE is a poly- 
chromatic Monte Carlo code that uses a 3-D adaptive 
grid to treat arbitrary geometries of emitting and ab- 
sorbing/scattering media. It self-consistently treats dust 
re-emission and self-absorption by iterating until the dust 
temperature converges. This is particularly important 
in the highly optically-thick central regions of merging 
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galaxies, and allows for future comparisons with longer 
wavelength observations by MlPS/Spitzcr, SCUBA (see, 
e.g., Narayanan et al. 2009), and the upcoming Her- 
schel mission. The latest version (Jonsson, P., Groves, 
B. & Cox, T. J. in prep.) includes a sophisticated pho- 
toionization code MAPPINGS and subrcsolution model 
to account for attenuation in the HII and photodissoci- 
ation regions around young stars (Groves et al. 2008). 
Outside these birth clouds, SUNRISE follows the light 
traveling through the hot phase of the ISM, assuming a 
negligible volume filling factor of the cold phase clouds, 
as was done in our line-of-sight attenuation code. In ad- 
dition to a different treatment of the radiative transfer 
(e.g., the inclusion of scattering), SUNRISE + MAP- 
PINGS also make use of a different stellar population 
synthesis code (Starburst99 by Leitherer et al. 1999). 

We ran SUNRISE on the simulations used in this pa- 
per, using a SMC dust model and identical initial con- 
ditions as adopted for our line-of-sight attenuation cal- 
culations. We adopted a birth cloud covering fraction 
of 0.3, but note that increasing the covering fraction to 
unity does not significantly alter our results. Populating 
our model universe using the merger rate function of §5, 
we obtain the model color distribution for the redshift 
interval 1.5 < z < 3 shown in Fig. 14. The color distri- 
bution is plotted separately for the runs with gas fraction 
/gas = 0.8 and / gas = 0.4 (see also §11.5), and includes 
both the pre- and the post-quasar phase. For reference, 
we also include panels illustrating the color distribution 
of massive observed galaxies in the same redshift inter- 
val (with masses derived from BC03 and M05 models), 
and realizations of the model photometry using our line- 
of-sight attenuation code in combination with BC03 and 
M05 stellar population synthesis templates. In each of 
the panels, the color-coding indicates the specific star for- 
mation rate. Typical error bars in the UVJ diagram of 
observed galaxies are drawn for the quiescent, and blue 
and red star-forming galaxies separately (their respec- 
tive regions in UVJ space are outlined by the polygons). 
Apart from the default SUNRISE photometry, we also 
present realizations in which one of the aspects of the 
radiative transfer is switched off. This allows us to dis- 
entangle the impact of scattering, birth clouds, and the 
multi-phase breakdown of the ISM. 

The default SUNRISE colors lie blueward of our line- 
of-sight attenuation photometry of the same set of sim- 
ulations, hence articulating the difficulty of reproducing 
the colors of observed dusty starbursts (objects in the 
dark-gray polygon of the observed panels in Fig. 14). 
Whereas -for the colors presented in this diagram- the 
birth cloud model by Groves et al. (2008) has a negligi- 
ble impact, we note that accounting for scattering causes 
a significant part of the blueward shift. Making the as- 
sumption that all gas is in the hot phase (i.e., switch- 
ing off the multi-phase breakdown) results in colors that 
are slightly redder during the star-forming phases than 
the default SUNRISE photometry, but still significantly 
bluer than the observed dusty starbursts. 

Our analysis highlights the need to understand the de- 
tails of translating physical quantities to observables (i.e., 
stellar population synthesis and radiative transfer) in or- 
der to effectively constrain models. Moreover, we note 
that -apart from the freedom in initial gas fraction / gas , 
which for our isolated mergers is set by hand- this issue 



applies equally, or perhaps more severely, to cosmologi- 
cal simulations, where the spatial resolution is often too 
limited to apply full radiative transfer. 

11.4. Dust distribution 

Possibly, the distribution of dust in the simulated 
galaxies might not reflect reality. A more efficient red- 
dening would be obtained if a foreground screen of ob- 
scuring material were present. One possible mechanism 
that could produce such a configuration on a galac- 
tic scale is a large-scale wind. The GADGET-2 code 
(Springcl 2005b) used to run the simulations in principle 
allows for such a phenomenon, but an investigation of 
the velocity field of the gas in the simulations is required 
to check whether such a wind is effectively taking place. 

On much smaller scales, below the resolution of the 
simulations, a more effective reddening might be ex- 
pected from taking into account that the molecular 
clouds in which new stars are formed, have a finite life- 
time. So far, we ignored the role of birth clouds (ex- 
cept for the comparison with SUNRISE radiative trans- 
fer, see §11.3), essentially assuming that they are instan- 
taneously dispersed when new stars form. Modeling the 
attenuation by HII and HI regions around young stars 
and by the ambient ISM, Chariot & Fall (2000) came up 
with a simple recipe where the attenuation of light from 
young stars is increased threefold with respect to the 
attenuation of the light from old stars, until the birth 
clouds disperse after ~ 10 7 yr. However, applying this 
recipe as a sub-grid model, the integrated SEDs of the 
simulated galaxies do not always redden over the entire 
optical-to-NIR wavelength range. 

We find that, until the star formation rate drops 
shortly after the final coalescence, a threefold extinc- 
tion toward young stars (< 10 Myr) reddens the inte- 
grated U — V color by - 0.15 mag. During the early 
merger phases, we observe a similar trend for the inte- 
grated V — J color, with birth clouds causing a redden- 
ing of - 0.1 mag. However, during the nuclear star- 
burst, the trend is reversed and integrated V — J colors 
bluen by - 0.15 mag. The reason is that, although 
young stars are intrinsically bluer, the column densities 
to the nuclear region where the starburst takes place 
are so much higher that, even ignoring the attenuation 
by birth clouds, the young component to the integrated 
light is redder in V — J than the old component. The 
balance between reddening and dimming by birth clouds 
subsequently reddens the already red young component, 
but downweights its contribution to the integrated light, 
thus producing a bluer overall color. At later times, the 
contribution of young stars is negligible, and so are the 
changes in the integrated colors when applying the sim- 
ple recipe for birth clouds. 

An alternative implementation, in which birth clouds 
around young stars (< 10 Myr) have a fixed Ay = 3 
instead of additional attenuation proportional to the line- 
of-sight column density through the entire galaxy, suffers 
from the same down-weighting of the young component 
to the integrated light. Simple birth cloud models of 
the kind described in this section seem insufficient to 
produce the colors of dusty starbursts. As described in 
§11.3, the more sophisticated birth cloud model adopted 
by SUNRISE causes an effective reddening of the U — V 
and V — J colors, but is also insufficient to account for 
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Fig. 14. — Rest-frame UVJ color-color diagram of massive galaxies at 1.5 < z < 3 as observed (big panels), and modeled with a variety 
of stellar population synthesis codes and radiative transfer methods (small panels). Median error bars are plotted for observed galaxies 
enclosed by the white, light-gray and dark-gray polygons. Massive galaxies exhibit a large range of specific star formation rates, both in 
the observations and simulations. The model colors depend on input stellar population synthesis, method of radiative transfer, and initial 
conditions of the simulations (e.g., initial gas fraction). The model colors of quiescent galaxies are generally consistent with those of their 
observed counterparts. However, in none of the realizations of the model photometry, the locus of star-forming galaxies extends as far 
redward (in U — V and particularly V — J) as the observations. 



the colors of observed dusty starbursts. 

It is very well possible that the detailed dust distri- 
bution, on scales far below the resolution of our sim- 
ulations, plays a crucial role in determining the galaxy 
colors. The processes that govern the distribution of dust 
are highly complex, involving wind escape from massive 
stars, dust grain formation, and diffusion in the interstel- 
lar medium. Unfortunately, we are not able to address 
those issues with self-consistent numerical modeling in 
our galactic-scale simulations, at the moment. Although 
the discussion in §11.3 and this subsection does not do 
justice to the full complexity of the problem, it serves to 
highlight the uncertainties in the simulated colors. 

11.5. Merger parameters 

In addition, the discrepancy in colors might imply that 
the simulations are not characteristic for the merger ac- 



tivity occuring in the real universe. E.g., from the shape 
of the stellar mass function it can be expected that galax- 
ies at the high-mass end are more likely to merge with 
galaxies of lower rather than comparable mass (Khochfar 
& Silk 2006; Peng 2007). Hopkins et al. (2006b) con- 
firmed the robustness of the model for quasar lifetimes 
and the derived merger rate function against changes in 
various parameters of the merging galaxies, such as gas 
fraction, orbital parameters and changes in the mass ra- 
tio of the progenitors (considering 1:1, 2:1, 3:1, and 5:1 
mass ratios). The conversion from a quasar birth rate 
to a spheroid birth rate relies on a proper knowledge of 
the black hole - bulge mass relation and how it evolves 
with redshift. This source of uncertainty will be reduced 
as more stringent observational constraints of the scaling 
relation at high redshift become available. Considering 
progenitor mass ratios, Dasyra et al. (2006) find for a 
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population of local ULIRGs that still have 2 distinct nu- 
clei that the typical mass fraction is 1.5:1, close to equal- 
mass mergers. In order to refine the model predictions, 
a detailed study of minor merger simulations is required 
to determine the minimum mass ratio required to trigger 
a (low-luminosity) quasar phase. 

Finally all of the simulations used in this work are 
equal-mass gas-rich mergers (/ gas = 0.8 at the start of 
the simulation). In order to study the dependence of our 
results on the adopted initial gas fraction / gas , we re- 
peated our analysis with / gas = 0.4 instead of / gas = 0.8. 
Although further gas infall is still not modeled, this test 
gives a crude but illustrative characterization of the de- 
pendence on the available gas reservoir and may repre- 
sent differences induced by the gas accretion or evolu- 
tionary history. The / gas = 0.8 runs represent a sce- 
nario in which most of the stars are formed during the 
merger. In the / gas = 0.4 runs on the other hand, the 
progenitor disks are more mature and a relatively smaller 
fraction of the final stellar mass is formed in episodes of 
merger-triggered star formation. The star formation his- 
tory prior to the start of the simulation is assumed to 
be constant, and the initial metallicities of gas and stars 
are set according to a closed box model (see Wuyts et al. 
2009 for more details). As a result, the / gas = 0.4 runs 
do not only start with a larger fraction of the baryonic 
mass in stars, but also with more metal-rich gas (0.9Z Q 
at the start of the simulation) and an initial stellar popu- 
lation having an older mass- weighted age (on the order of 
1 Gyr) and larger mass-weighted metallicity (~ 0AZ Q ). 

In broad terms, our conclusions remain unaltered when 
lowering the adopted / gas from 0.8 to 0.4. Especially 
the results regarding the post-quasar phase seem robust, 
since by that time the stellar mass evolution of the low 
and high gas fraction runs have converged and most of 
the sensitivity to the initial conditions has been washed 
out. In other words, the post-quasar systems have dis- 
tinctly lower SFR/M values than the earlier evolution- 
ary stages, their UVJ colors lie in the selection wedge 
that was designed to select the observed quiescent galax- 
ies, and, as for the / gas = 0.8 results, we find that the 
model prediction for their number and mass density as 
derived from the observed quasar luminosity function ex- 
ceeds the observed abundance by a factor of ~ 2.5. 

For the model predictions of the star-forming popu- 
lation, the following dependencies on / gas are observed. 
First, in the / gas = 0.4 runs more stars already formed 
prior to the merger-triggered star formation phase. Some 
systems in the early stages of merging, that in the 
/ g as = 0.8 realization had not built up a sufficient amount 
of mass to enter the sample, will now fall above the mass 
limit, hence boosting the model prediction and reducing 
the discrepancy with respect to the abundance of ob- 
served massive star-forming galaxies by a factor ~ 1.7. 
Second, the initial conditions affect the intrinsic photom- 
etry of simulated merging galaxies by the increased age 
and metallicity of the stellar population, and the atten- 
uated photometry is affected also by the smaller amount 
and larger metal-content of the gas. The presence of an 
old underlying population and relatively smaller build-up 
of new stars during the merger explains why the tail to- 
wards blue V— J and especially blue U—V colors (see Fig. 
6 and Fig. 7) is lacking in the color distribution realized 
with / gas = 0.4. For our default BC03 models, we find 



the U — V color distribution for the entire (pre- and post- 
quasar) population to range from 1 to 2 with the central 
68% interval bracketed by 1.3 < U - V < 1.8. The V — J 
distribution ranges from 0.5 to 1.6, with a narrow central 
68% interval of 0.9 < V - J < 1.1. The dependence on 
population synthesis code and attenuation law is similar 
as for the higher / gas runs. E.g., using M05 models we 
find a V — J distribution ranging from 1 to 2 (central 
68% interval 1.3 < V - J < 1.4). We conclude that the 
discrepancy with respect to the optical-to-NIR colors of 
observed dusty starbursts is not alleviated by adopting 
a different initial gas fraction. This is also illustrated in 
Fig. 14. 

11.6. Evolutionary history 

Alternatively, it is possible that dusty starburst galax- 
ies are not triggered by mergers, but had a different evo- 
lutionary history. Daddi et al. (2007a) make this claim 
based on the long star formation timescales of ULIRGs 
at high rcdshift, and the relatively tight relation between 
SFR and stellar mass. We note, however, that our sim- 
ulations of isolated disk galaxies with initial conditions 
identical to those of the merger progenitors also fail to 
produce actively star-forming systems with colors simi- 
lar to dusty red starbursts. While CO interferometry by 
Tacconi et al. (2006, 2008) shows key evidence that ma- 
jor merging is taking place in essentially all SMGs, the 
merger fraction among the general z ~ 2 star-forming 
population may be lower. Spatially-resolved kinematic 
studies by the SINS survey (Forster Schreiber et al. 
2006b; Genzel et al. 2006, 2008) and van Starkenburg 
et al. (2008) show a significant number of rotating disks 
with high (~ 100M Q yr _1 ) SFR but no evidence for re- 
cent or on-going major merging. Applying the method of 
kinemetry to the SINS sample, Shapiro et al. (2008) finds 
a third of z ~ 2 star-forming galaxies to be undergoing 
major merging, consistent with our model prediction. 

From a theoretical perspective, recent cosmological hy- 
drodynamic simulations have suggested that most of the 
gas accretion at the high-mass end takes place in cold 
flows along dark matter filaments, with the streams being 
~ 50% smoothly flowing material, and the other ~ 50% 
in clumps of mass ratio < 10 : 1 with respect to the 
accreting system (Keres et al. 2005; Dckel & Birnboim 
2006; Ocvirk et al. 2008; Dekel et al. 2009). Unlike 
major mergers, such flows could keep the rotating disk 
configuration observed in many star-forming z ~ 2 galax- 
ies intact. Using a semi- analytic model with star forma- 
tion and feedback recipes based on hydrodynamic simu- 
lations, Somerville et al. (2008) also find that most of the 
global star formation occurs in a quiescent mode, rather 
than in merger-induced starbursts. This could explain 
why our model prediction for the abundance of merger- 
triggered star-forming galaxies accounts for only a third 
of all observed star-forming galaxies at the high-mass end 
(see §9.2 and §11.1.2). 

We conclude that, while our results are consistent 
with a scenario where all massive quiescent galaxies at 
1.5 < z < 3 have formed by a merging event that trig- 
gered quasar activity, it leaves room for other mecha- 
nisms than major mergers contributing significantly to 
the z ~ 2 star-forming population. It would be interest- 
ing to scrutinize such alternative scenarios in the same 
way as presented here. 
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11.7. Mass loss and intergalactic environment 

Gas replenishment from mass loss and infall of gas from 
the intergalactic environment could change the optical 
depths and thus the reddening factors. The simulations 
only take into account a small amount of mass loss: 10% 
of the gas mass converted into stars is instantaneously re- 
turned to the interstellar medium, accounting for short- 
lived stars that die as supernovae (Springel & Hernquist 
2003). The total fraction of the mass lost by an aging sin- 
gle stellar population with Salpcter (1955) IMF amounts 
to <~ 30% (BC03) and is even higher (~ 50%) for more re- 
alistic IMFs such as Kroupa (2001) or Chabrier (2003). 
Furthermore, the simulations do not allow for infall of 
primordial gas at later times. Consequently, they can- 
not prove that descendants of galaxies that once showed 
up in the quasar luminosity function and after the shut- 
down of star formation reached red colors, will remain 
quiescent forever. Small amounts of newly accreted gas 
triggering star formation may be enough to shift a post- 
quasar galaxy outside the quiescent region of color-color 
space defined by Eq. 5, thus dropping their contribution 
to the observed galaxy population of massive quiescent 
red galaxies. Cosmological simulations at sufficient res- 
olution might resolve this problem. At the very least, 
it would be interesting to test the behavior of simulated 
merger remnants hosting a supermassive black hole when 
a small but continuous gas supply is applied. 

We note that a rejuvenation of part of the post-quasar 
population would simultaneously improve the agreement 
with the observations for both the star-forming and the 
quiescent galaxy abundance. 

11.8. Cosmic variance 

From the observational side, cosmic variance is the 
dominant source random of uncertainties for the determi- 
nation of the number and mass density of massive galax- 
ies. At z < 2.5, this is even the case when including the 
wider area MUSYC survey (Marchesini et al. 2008). This 
is particularly true for the quiescent population, that 
shows a stronger clustering than the star-forming one 
(Williams et al. 2009) . Surveys over a significantly larger 
area than FIRES and GOODS-South, but probing simi- 
lar depths at optical-to-MIR wavelengths, arc required to 
better constrain the fraction of massive quiescent galax- 
ies that post-quasar galaxies can account for. The ultra- 
VISTA survey by Dunlop, Franx, Fynbo,& LeFevre will 
provide NIR imaging over half of the 2 square degrees 
COSMOS field to an unprecedented depth. In combi- 
nation with very deep IRAC imaging during the warm 
Spitzer mission, such multi-wavelength surveys will allow 
to simultaneously probe further down the mass function 
of quiescent galaxies, and reduce the effect of cosmic vari- 
ance. 

12. SUMMARY 

We confronted the model by Hopkins et al. (2006b) 
with observations of massive galaxies at 1.5 < z < 3. The 
model translates the observed quasar luminosity func- 
tion into the abundance of massive merging galaxies and 
merger remnants. We derived the synthetic photome- 
try for these systems from a set of binary merger SPH 
simulations by Robertson et al. (2006a, 2006b) and T. 
J. Cox, with a range of masses, and including stellar 



and AGN feedback. We extracted mass-limited sam- 
ples of 1.5 < z < 3 galaxies with M > 4 x 10 10 M Q 
and M > 10 11 Mq from the FIRES+FIREWORKS 
and FIRES+FIREWORKS+MUSYC surveys respec- 
tively. We tested the model by comparing the predicted 
number and mass densities, the U — V and V — J color 
distributions, and SFR/M versus mass distribution with 
our observations of massive galaxies at 1.5 < z < 3. 

We find that the overall number density of galaxies 
with M > 4 x 10 10 M Q in the FIRES and GOODS-South 
fields (n = i.OtH x 10_4 Mpc~ 3 at z ~ 1.9 and n = 
3.3^2 x 10~ 4 Mpcr 3 at z ~ 2.6) is consistent within the 
uncertainties with the model prediction (n = 4.8 — 5.5 x 
10~ 4 Mpc~ 3 at z ~ 1.9 and n = 1.9 - 2.2 x 10~ 4 Mpc~ 3 
at z <~ 2.6). Likewise, the results obtained for the mass 
density are consistent: p* = 3.8l 1 ' 2 ! x 10 7 Mq Mpc~ 3 at 
z ~ 1.9 and p* = x 10 7 M Mpc~ 3 at z ~ 2.6 for 

the observations and = 5.4 — 6.0 x 10 7 M Q Mpc~ 3 at 
z ~ 1.9 and p* = 1.9 - 2.2 x 10 7 M Mpc~ 3 at z ~ 2.6. 

Separating massive galaxies by type, we find that the 
model photometry of the post-quasar population coin- 
cides with the region of U — V versus V — J color-color 
space that was defined by Labbe et al. (in prepara- 
tion) to select quiescent red galaxies. The modeled num- 
ber and mass densities of massive (M > 4 x 10 10 Mq) 
quiescent galaxies is consistent with the observations at 
z ~ 2.6, but somewhat larger (2-3 times) than observed 
at z <~ 1.9. The results based on the UVJ diagnostic 
diagram and on a MlPS-based SFR/M threshold are 
consistent. 

We added the MUSYC survey to our sample, increas- 
ing the area by a factor of 3.6, but by the shallower depth 
restricting our analysis to M > 10 11 Mq galaxies. From 
this sample, we derive qualitatively similar results. 

Although less constrained, the predicted abundances of 
galaxies with merger-triggered star formation (according 
to their UVJ colors or SFR/M > 1/iHubbic can also ac- 
count for a significant fraction of the observed actively 
star-forming galaxies (one third when using masses based 
on BC03 templates and taking into account biases in 
SED modeling). However, the predicted color distribu- 
tion of star-forming galaxies does not match the observa- 
tions. In particular the colors of red (V — J > 1.8) dusty 
starburst galaxies are not reproduced. We suggest a 
number of explanations for the lack of dusty red starburst 
galaxies in the model predictions. Possible reasons are 
an incomplete simulation of the observing procedure, dif- 
ferences in stellar population properties or merger char- 
acteristics between the observed and simulated galaxies, 
a different history for dusty starbursts than a merger- 
triggered scenario, infall of additional gas (and dust) 
from the intergalactic environment or mass loss, and a 
different distribution of the dust, e.g., caused by the pres- 
ence of large-scale outflows or birthclouds around young 
stars. 

Finally, we find hints of a pair excess at small angular 
scales, further strengthening the hypothesis that mergers 
play a key role in galaxy evolution. 

We conclude that the star formation in remnants of 
merger simulations is quenched abruptly, leading to col- 
ors that correspond well to those of observed massive 
galaxies with a quiescent stellar population. Using a 
merger rate derived from the observed quasar luminos- 
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ity function, we obtain number and mass densities of the 
quiescent population at 1.5 < z < 3 that are consistent 
within a factor of 2 to 3 with the observations. Possi- 
bly, the overprediction at z ~ 1.9 suggests the need to 
include gas infall as refinement to the model. The pre- 
dicted abundance of merger-triggered star-forming sys- 
tems accounts for 30-50% of the observed star-forming 
population, leaving ample room for other star formation 
mechanisms than major merging. The most serious chal- 
lenge to the model is posed by the color distribution of 
star-forming galaxies, which is not well reproduced. The 



detailed dust distribution, on a galaxy- wide scale and/or 
scales far below the resolution of the simulation, may well 
cause this problem. With this work, we hope to encour- 
age further investigations focussing on the translation of 
simulated physical quantities to observables. 
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